LIVESENSE Data Analytics Blog

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

LivesenseDataAnalyticsBlog

MENU

Factorization Machinesをレコメンデーションで使うときの評価推定値計算

こんにちは、リブセンスで統計や機械学習関係の仕事をしている北原です。今回はレコメンデーションで使う評価推定値計算の効率化に関する小ネタです。機械学習を実務で使うときのちょっとした工夫に関するお話です。実装にはJuliaを使います。

FM(Factorization Machines)をレコメンデーションで使う場合、各ユーザーに対してレコメンド可能なアイテムの評価推定値計算を行うため、ユーザー数とアイテム数が多くなると非常に計算時間がかかります。学習データで交差検証をしたりコンペで使ったりしているだけだとあまり問題にならないのですが、実務では結構問題になります。こういうのは実務の現場で個別に対処されているためか知見もあまり公開されていないようです。対処方法はいろいろあるのですが、ここでは計算方法を少し工夫することで計算量を削減する簡単な方法を紹介します。合わせて計算量の見積もりや計算量削減の考え方、アルゴリズム実装のコツについてもお伝えできればと思います。

ユーザーベース協調フィルタリングの評価推定値計算

ユーザーベース協調フィルタリングを利用してレコメンデーションを行う場合、ユーザーとアイテムの組み合わせごとに評価推定値を計算し、各ユーザーに対して評価推定値が高いアイテムをレコメンドするのが一般的です。そのため、ユーザー数やアイテム数が多いと計算量が多くなり、実務でレコメンデーションを利用するときに計算時間が問題になることがよくあります。交差検証で精度が高そうだったから運用にのせようと考えて、一通りのレコメンデーション処理をやってみたら推定値計算に時間がかかりすぎることに気がついた、といったこともあるのではないでしょうか。学習は高速なので交差検証などを行っている範囲では計算時間の問題に気づきにくいのです。

まず、評価推定値計算に必要なデータ量について考えてみましょう。弊社サービスの場合、ユーザー数、アイテム数とも、多くとも数十万程度を見込んでおけばよいのでデータ量はそれぞれ$O(10^5)$です。したがって、組み合わせの数は$O(10^{10})$となります。つまり、評価推定値計算のデータ量は$O(10^{10})$です。評価データとして、口コミ評点や応募、求人詳細ページ閲覧を使った場合、1ユーザーあたりの評価データ数の平均はせいぜい$O(10)$ぐらいです。そのため、学習データとして使われる評価データ数のオーダーは$O(10^6)$になります。学習時と推定値計算時とでは1サンプルデータあたりの計算量が異なるので単純比較はできないものの、データ数が一万倍も多いのであれば、学習時だけでなく推定値計算時の計算量についても配慮が必要ではないかと予想がつきます。レコメンデーションを使いたい場面では、多くの場合1ユーザーあたりの評価データ数はユニークアイテム数より十分小さいので、評価推定値の計算量は弊社に限らず問題になりやすいのではと思います。

対処方法はいろいろあるものの、代表的な方法は、運用の工夫、並列・分散処理の利用、計算量改善の三つです。どれか一つしか使えないわけではなく組み合わせて使うこともできるので、組み合わせるのが一般的です。

運用の工夫としては、レコメンデーション効果の低い組み合わせを計算しないようにする方法があります。例えば、求人サイトの場合、会員登録時に希望職種や希望勤務地を入力してもらうことが多いので、希望と合致する求人しかレコメンドしないように制限する方法が考えられます。この方法は簡単でよく使われているのですが、希望職種や希望勤務地が適切に登録されてないと見逃しが発生してしまうという問題があります。同じような業務内容なのに異なる職種として登録されている求人はレコメンドされませんし、希望勤務地の隣接地域に非常にマッチする求人があってもレコメンドされません。制限を強めると見逃しが多くなる一方で、弱めると組み合わせが多くなるのでバランスをとる必要があります。

並列・分散処理の利用はいうまでもないですね。単純に組み合わせが多いだけで独立して計算可能な部分が多いので並列・分散処理を使うのは容易です。高性能なインスタンスや分散処理システムを利用するなどスケールさせるのも比較的簡単です。しかし、並列・分散処理だけで解決しようとすると計算機リソースを消費しすぎてコストがかかるという問題もあります。

計算量改善は計算式を整理して計算量を削減する方法です。計算式をうまく整理できると実装変更以外のコストをかけずに計算時間を短縮することができます。ただし、アルゴリズムやモデルによっては使えないこともあります。シミュレーションなどの数値計算や、最適化、機械学習などのアルゴリズム実装をすることが多い人にとっては、計算式を整理して効率的な計算ができるようにするのは当たり前のことなので、ノウハウが公開されることは少ないです。

Factorization Machinesの評価推定値計算

それでは、FMで評価推定値計算するときの計算量について考えてみましょう。

まず、評価推定値計算時のFMについて説明します。学習時は評価済みデータからモデルパラメータを推定しますが、評価推定値計算時はユーザーコンテキストとアイテムコンテキストを組み合わせたデータの評価推定値を計算します。ユーザー数を$s^u$、アイテム数を$s^i$とすると、$s=s^u s^i$の評価推定値を計算します。全コンテキスト数(特徴量数)を$n$、因子数を$k$とし、評価推定値$\hat{\mathbf{y}} \in {\mathbb R}^{s}$、評価と対応するコンテキストデータ${\bf X} \in{\mathbb R}^{s \times n}$(1サンプルのコンテキストは$\bf{x}$と表記)、モデルパラメータ$w_0 \in {\mathbb R}$、${\bf w} \in {\mathbb R}^n$、${\bf V} \in {\mathbb R}^{n \times k}$を使って、評価推定値$\hat{y}$は

$$ \begin{eqnarray} \hat{y}({\bf x}) &=& w_0 + \sum^n_{j=1}{w_j x_j} + \sum^n_{l=1}\sum^n_{j=l+1}{\hat{w}_{l, j} x_l x_j} \label{eq:fm} \\ \hat{w}_{l, j} &=& \sum^k_{f=1}{v_{l,f} v_{j,f}} \label{eq:w} \end{eqnarray} $$

と計算できます。評価推定値計算時にはモデルパラメータ$w_0$、$\bf{w}$、$\bf{V}$は既知です。コンテキスト$\bf{X}$はユーザーに関するコンテキスト${\bf X}^u \in{\mathbb R}^{s \times n^u}$とアイテムに関するコンテキスト${\bf X}^i \in{\mathbb R}^{s \times n^i}$が連結されたもの${\bf X} = [{\bf X}^u \, {\bf X}^i]$になっています。ここで$n^u$と$n^i$はそれぞれユーザーとアイテムのコンテキスト数で$n=n^u+n^i$です。それに対応し、モデルパラメータも${\bf w}^u \in {\mathbb R}^{n^u}$、${\bf w}^i \in {\mathbb R}^{n^i}$、${\bf w} = [{\bf w}^u \, {\bf w}^i]$、${\bf V}^u \in {\mathbb R}^{n^u \times k}$、${\bf V}^i \in {\mathbb R}^{n^i \times k}$、${\bf V}=\begin{bmatrix}{\bf V}^u \ {\bf V}^i \end{bmatrix}$と表記します。

このまま計算すると計算量は$O(s^u \, s^i \, k \, (n^u + n^i)^2)$になりますが、交互作用項を変形した

$$ \begin{eqnarray} \hat{y}({\bf x}) &=& w_0 + \sum^n_{j=1}{w_j x_j} + \frac{1}{2} \sum_{f=1}^{k}\left(\left(\sum_{j=1}^{n} v_{j, f} x_{j}\right)^{2}-\sum_{j=1}^{n} v_{j, f}^{2} x_{j}^{2}\right) \label{eq:fm2} \end{eqnarray} $$

を使うと、計算量は$O(s^u \, s^i \, k \, (n^u + n^i))$になります。FMでレコメンデーションを行うときのコンテキストデータはスパースなので、1ユーザーあるいは1アイテムあたりのコンテキストの非ゼロデータ個数を$\overline{m}^u$、$\overline{m}^i$とすると、計算量は$O(s^u \, s^i \, k \, (\overline{m}^u + \overline{m}^i))$となります。

では具体的な計算量がどのくらいになるか考えてみましょう。弊社の場合だと、ユーザー数とアイテム数はそれぞれ数十万程度を見込んでおけばよいので$s^u \sim O(10^5)$、$s^i \sim O(10^5)$、因子数は多くとも数百程度とすると$k \sim O(10^2)$、コンテキストの非ゼロデータ個数は、会員登録情報など簡単な前処理だけで使えるデータであればせいぜい数十、求人詳細テキストなどを使う場合でも数百程度に収まると考えると$\overline{m}^u + \overline{m}^i \sim O(10^2)$になります。そうすると、計算量は$O(s^u \, s^i \, k \, (\overline{m}^u + \overline{m}^i)) \sim O(10^{14})$になります。かなり大きそうですが比較対象がないとわかりにくいですね。

そこで、学習時と比較してどのぐらいの計算量になるか考えてみましょう。以前の記事でFMのモデルパラメータをAlternating Least Squares(ALS)で推定する方法を紹介したので、それと比較します。学習時の計算量は、評価数$s$、因子数$k$、1モデルパラメータに対応する非ゼロデータ個数の平均を$\overline{m}_C$とすると$O(s \, k \, \overline{m}_C)$でした。弊社サービスの場合、1ユーザーあたりの平均評価数は$O(10)$程度なので$s \sim O(10^6)$ になります。ユーザーIDやアイテムID以外のコンテキストが極端に多くなければコンテキスト数$n$はユーザー数$s^u$やアイテム数$s^i$とほぼ同じオーダーになるので$n \sim O(10^5)$です。1モデルパラメータに対応する非ゼロデータ個数$\overline{m}_C$と1データあたりのコンテキストの非ゼロデータ個数$\overline{m}^u + \overline{m}^i$との比は、評価数の逆数$1/s$とコンテキスト数の逆数$1/(n^u + n^i)$との比と等しくなるので、$\overline{m}_C \sim O(10)$になります。そのため、学習時の計算量は$O(s \, k \, \overline{m}_C) \sim O(10^{9})$となります。つまり、弊社の場合、学習時と比較して評価推定の計算量は十万倍程度になります。学習は一分で終わっても評価推定計算にはおよそ七十日かかるので、並列・分散処理で十倍程度高速化しても一週間はかかります。求人サービスの場合は求人の入れ替わりが激しいので、時間をかけて高精度なモデルを作るより頻繁にモデルを更新して新しい求人をレコメンドできるようにしたほうがよいことが多いです。そのため、これぐらい計算時間がかかると運用上のボトルネックになってしまいます。

特徴量作成を行う立場からも考えてみましょう。レコメンデーションの精度を上げるためには、有用なコンテキスト情報を増やす必要があります。求人サービスの場合、希望勤務地などの会員情報や年収などの求人情報だけでなく、求人詳細テキストから自然言語処理を利用して作成した特徴量やユーザーが過去に閲覧した求人や閲覧傾向の変化などレコメンデーションに有用と思われる情報は数多くあります。そのため、コンテキストを増やしていきたいですが、評価推定値の計算時間も増えてしまうので容易には増やせないというジレンマが生まれてしまいます。

事前計算による計算量削減

それでは、本記事のテーマである計算量削減方法について説明します。

基本的な考え方は非常にシンプルで、繰り返し同じ計算をする部分については一回計算したらその値を記憶しておき繰り返し使い回すことで計算量を削減します。式($\ref{eq:fm}$)に含まれる計算の中で繰り返し同じ計算をする部分というのは、ユーザーあるいはアイテムに関するモデルパラメータとコンテキストデータのみで計算できる部分です。

計算量削減のコツは課題を定式化して改善可能な部分を明確にすることです。そこで、まずユーザーとアイテムの違いがわかるように式($\ref{eq:fm}$)を書き直してみましょう。前節で導入した表記を使うと式($\ref{eq:fm}$)は

$$ \begin{eqnarray} \hat{y}({\bf x}) &=& w_0 + \sum^{n^u}_{j=1}{w^u_j x^u_j} + \sum^{n^i}_{j=1}{w^i_j x^i_j} + \sum^{n^u}_{l=1}\sum^{n^u}_{j=l+1}{\hat{w}^{u}_{l, j} x^u_l x^u_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:fm3} \\ \hat{w}^u_{l, j} &=& \sum^k_{f=1}{v^u_{l,f} v^u_{j,f}} \label{eq:w_u} \\ \hat{w}^i_{l, j} &=& \sum^k_{f=1}{v^i_{l,f} v^i_{j,f}} \label{eq:w_i} \\ \hat{w}^{u,i}_{l, j} &=& \sum^k_{f=1}{v^u_{l,f} v^i_{j,f}} \label{eq:w_ui} \end{eqnarray} $$

と書くことができます。和の記号の範囲や上付き添字で示されたユーザーとアイテムの違いに注意してください。式($\ref{eq:fm3}$)の右辺第6項のみがユーザーとアイテムの組み合わせになっており、それ以外はユーザーのみあるいはアイテムのみのデータで計算可能で、事前に一回計算して結果を記憶しておけば$O(1)$で計算できることがわかります。式($\ref{eq:fm3}$)の右辺第6項をそのまま計算すると計算量は$O(k \, n^u \, n^i)$、スパースなことを考慮すると$O(k \, \overline{m}^u \, \overline{m}^i)$となり計算量が小さいとは言い難いので、ここが改善ポイントではないかと予想がつきます。

念のため、式($\ref{eq:fm3}$)の第4項、第5項の計算量についても考えておきましょう。式($\ref{eq:fm3}$)の第4項、第5項は、式($\ref{eq:fm}$)の交互作用項と同じ形をしているので同様の変形をすると、式($\ref{eq:fm2}$)の交互作用項と同じ形になります。式($\ref{eq:fm3}$)の第4項は各ユーザーについて、第5項は各アイテムについて計算する必要があるので、計算量はそれぞれ$O(s^u \, k \, \overline{m}^u)$と$O(s^i \, k \, \overline{m}^i)$になります。前節と同様に$s^u \sim O(10^5)$、$k \sim O(10^2)$、$\overline{m}^u \sim O(10^2)$とすると、式($\ref{eq:fm3}$)の第4項の計算量は$O(s^u \, k \, \overline{m}^u) \sim O(10^9)$となるため学習時より計算量が小さいことがわかります。式($\ref{eq:fm3}$)の第5項も同様なので、学習時に計算量が問題にならなければ、式($\ref{eq:fm3}$)の第4項、第5項でも問題になることはないと考えてよいでしょう。

次に、式($\ref{eq:fm3}$)の第6項について考えましょう。式($\ref{eq:fm3}$)の第6項は

$$ \begin{eqnarray} \sum^{n^u}_{l=1}\sum^{n^i}_{j=1}{\hat{w}^{u,i}_{l, j} x^u_l x^i_j} &=& \sum^{n^u}_{l=1}\sum^{n^i}_{j=1}\sum^k_{f=1}{v^u_{l,f} v^i_{j,f} x^u_l x^i_j} \nonumber \\ &=& \sum^k_{f=1} \Biggl(\sum^{n^u}_{l=1} {v^u_{l,f} x^u_l} \Biggr) \Biggl(\sum^{n^i}_{j=1} {v^i_{j,f} x^i_j} \Biggr) \nonumber \\ &=& \sum^k_{f=1} q^u_f q^i_f \label{eq:qq} \end{eqnarray} $$

と変形することができます。ここで

$$ \begin{eqnarray} q^u_f &=& \sum^{n^u}_{j=1} {v^u_{j,f} x^u_j} \label{eq:q_u} \\ q^i_f &=& \sum^{n^i}_{j=1} {v^i_{j,f} x^i_j} \label{eq:q_i} \end{eqnarray} $$

としています。$q^u_f$はユーザーデータのみで計算でき、$q^i_f$はアイテムデータのみで計算できるので、事前に一回計算し記憶しておけば式($\ref{eq:qq}$)の計算量は$O(k)$になります。したがって、式($\ref{eq:fm3}$)の計算量は$O(s^u \, s^i \, k)$となります。

もう少し$q^u_f$と$q^i_f$について考えてみましょう。$q^u_f$はユーザーごとに因子数分の$k$個計算しておく必要があります。つまり、合計$k \times s^u$個計算しておく必要があるので、因子を行、ユーザーを列とした$q^u_f$を要素にもつ${\bf Q}^u \in {\mathbb R}^{k \times s^u}$を事前に計算し記憶しておきます。同様に、アイテムについても因子を行、アイテムを列とした$q^i_f$を要素にもつ${\bf Q}^i \in {\mathbb R}^{k \times s^i}$を計算し記憶しておきます。${\bf Q}^u$を作成する計算量は$O(s^u \, k \, \overline{m}^u)$、${\bf Q}^i$を作成する計算量は$O(s^i \, k \, \overline{m}^i)$で、式($\ref{eq:fm3}$)の第4項、第5項と同じ計算量になるため、計算量が問題になることはありません。

まとめると、ユーザーのみあるいはアイテムのみで計算可能な部分を事前に計算しておき、その結果を式($\ref{eq:fm3}$)の計算時に使うことで式($\ref{eq:fm3}$)の計算量は$O(s^u \, s^i \, k)$となります。それにより、コンテキストの非ゼロデータ個数$\overline{m}^u$や$\overline{m}^i$には依存せず、コンテキストを増やしても評価推定値計算の計算量は増加しなくなります。前節と同様の想定で$s^u \sim O(10^5)$、$s^i \sim O(10^5)$、$k \sim O(10^2)$とすると$O(s^u \, s^i \, k) \sim O(10^{12})$となります。並列・分散処理でさらに十倍程度高速化できれば計算量は$O(10^{11})$になるので、前節で見積もった学習時の計算量$O(10^9$)より百倍大きいものの、これぐらいであれば現実的な時間内に計算を終えることができるようになります。弊社でもここで紹介した計算方法と並列処理を併用して運用しています。

なお、式($\ref{eq:fm3}$)の計算量$O(s^u \, s^i \, k)$は、Matrix Factorization(MF)の評価推定値計算と同じ計算量になっています。そのため、レコメンデーションでMFを使った運用ができているならば、コンテキストを使うFMを導入したときの評価推定値計算時間は同程度になるため、計算時間がボトルネックになってFMを導入できなくなることはほとんどありません。

実装

では実装してみましょう。今回も簡潔な表記が可能なJuliaを使います。今回はコードも少ないので少し実装のコツも紹介します。

今回扱っている数式はシンプルなので実装するのは比較的容易だと思いますが、数式はわかるしコードも読めるけど慣れていないと自分で書くとうまく書けないという人もいるかもしれません。アルゴリズムの種類によって数式を実装に落とすコツは異なるのですが、ここでは多くの場合に使える基本的な考え方を2つ紹介します。

1つ目は、ベクトルや行列などの変数を作成するときのサイズと型についてです。論文などを読むと記号の定義が書いてあるのでそれを使います。例えば、${\bf V} \in {\mathbb R}^{n \times k}$となっていれば、行数$n$、列数$k$の浮動小数点数型の行列にすればよいことがわかります。そのため、記号を定義しているところを全て見つけてベクトルや行列を作成するコードを書いていけば、実装しようとしているアルゴリズムに必要な変数は一通りコードに落とせます。

2つ目は、和の記号についてです。和の記号を使った計算は、ループに置き換えることができるというのは基本ですよね。和の記号が三つかかっている計算は三重ループになります。注意するのは添字と和の範囲です。論文によっては和の範囲が省略されていたり実装しにくかったりするので、アルゴリズムを実装する前に確認しておく必要があります。例えば、$\sum^k_{f=1}$とあればループ変数$f$を1から$k$まで変えるループであるとすぐにわかりますが、$\sum_{f \in F}$とあれば集合$F$の定義を確認しておく必要があります。$F$の要素は特別な条件を満たさなければいけないのであれば条件処理の追加が必要になったりします。また、変数の添字が省略されていることもよくあります。論文では誤解されない範囲で簡潔な表記になっていて、実装には適していない表記になっていることもあるので、そのような場合は和の範囲を明示したり添字を省略しない形に書き直すと実装しやすくなります。

1データ分の推定値計算と、${\bf Q}^u$や${\bf Q}^i$の事前計算のコードは以下のとおりです。

using Distributions
using SparseArrays
using LinearAlgebra

# 以前の記事(https://analytics.livesense.co.jp/entry/2019/05/28/110000)で紹介したコードと同じ
function y_hat(x_, w₀, w, v)
    (i, x) = findnz(x_)  # 非ゼロデータのインデックスと値
    n = nnz(x_)          # 非ゼロデータの数
    k = size(v)[2]

    ΣΣwᵢⱼxᵢxⱼ = 0.0
    for f_ in 1:k
        Σvx = 0.0
        Σv²x² = 0.0
        for j_ in 1:n
            vx = v[i[j_], f_] * x[j_]
            Σvx += vx
            Σv²x² += vx^2
        end
        ΣΣwᵢⱼxᵢxⱼ += Σvx^2 - Σv²x²
    end
    return w₀ + x ⋅ w[i] + 0.5 * ΣΣwᵢⱼxᵢxⱼ  # 式(3)
end

function precompute_q(x_, v, s, k)
    q = zeros(Float64, k, s)  # (因子数, データ数) 

    for i_ in 1:s
        (idx, x) = findnz(x_[i_, :])
        nnz_ = nnz(x_[i_, :])

        for f_ in 1:k
            for j_ in 1:nnz_
                q[f_, i_] += v[idx[j_], f_] * x[j_]  # 式(9), (10)
            end
        end
    end
    return q
end

呼び出しコードは以下のとおりです。事前計算を使うと非ゼロデータ個数が増えても計算量がほとんど変わらないことを確認してみてください。

# サンプルデータ作成
sparseness = 0.01  # 非ゼロの比率

sᵘ = 10^3           # ユーザー数
nᵘ = 10^3           # ユーザーコンテキスト数
xᵘ = convert(SparseMatrixCSC{Int64,Int64}, sprand(Bool, sᵘ, nᵘ, sparseness))
xᵘᵗ = sparse(xᵘ')

sⁱ = 10^3           # アイテム数
nⁱ = 10^3           # アイテムコンテキスト数
xⁱ = convert(SparseMatrixCSC{Int64,Int64}, sprand(Bool, sⁱ, nⁱ, sparseness))
xⁱᵗ = sparse(xⁱ')


# モデルパラメータはすでに計算済みという想定(ここでは乱数を使う)
k = 10
n = nᵘ + nⁱ
w₀ = rand()
w = rand(n)
v = rand(Normal(0.0, 0.001), n, k)

# 事前計算用にユーザーとアイテムのパラメータを分ける
wᵘ = w[1:nᵘ]
vᵘ = v[1:nᵘ, :]
wⁱ = w[(nᵘ+1):n]
vⁱ = v[(nᵘ+1):n, :]


# 事前計算なし
@time e1 = [y_hat(vcat(xᵘᵗ[:,i], xⁱᵗ[:,j]), w₀, w, v) for j in 1:sⁱ, i in 1:sᵘ]

# 事前計算あり
@time ŷᵘ = [y_hat(xᵘᵗ[:, i], 0, wᵘ, vᵘ) for i in 1:sᵘ]                              # 式(4)の第2項、第4項
@time qᵘ = precompute_q(xᵘ, vᵘ, sᵘ, k)                                            # 式(9)
@time ŷⁱ = [y_hat(xⁱᵗ[:, i], 0, wⁱ, vⁱ) for i in 1:sⁱ]                                  # 式(4)の第3項、第5項
@time qⁱ = precompute_q(xⁱ, vⁱ, sⁱ, k)                                               # 式(10)
@time e2 = [w₀ + ŷᵘ[i] + ŷⁱ[j] + qᵘ[:,i] ⋅ qⁱ[:,j] for j in 1:sⁱ, i in 1:sᵘ]  # 式(4)、式(8)


# 事前計算が正しく動いているか確認
(e1 .- e2).^2 |> mean |> sqrt



# 非ゼロ比率を一定のままコンテキストデータサイズを増やす => 非ゼロデータ個数が増える
nᵘ = 10^4           # ユーザーコンテキスト数
xᵘ = convert(SparseMatrixCSC{Int64,Int64}, sprand(Bool, sᵘ, nᵘ, sparseness))
xᵘᵗ = sparse(xᵘ')

nⁱ = 10^4           # アイテムコンテキスト数
xⁱ = convert(SparseMatrixCSC{Int64,Int64}, sprand(Bool, sⁱ, nⁱ, sparseness))
xⁱᵗ = sparse(xⁱ')

n = nᵘ + nⁱ
w₀ = rand()
w = rand(n)
v = rand(Normal(0.0, 0.001), n, k)

wᵘ = w[1:nᵘ]
vᵘ = v[1:nᵘ, :]
wⁱ = w[(nᵘ+1):n]
vⁱ = v[(nᵘ+1):n, :]

# 事前計算なし
@time e1 = [y_hat(vcat(xᵘᵗ[:,i], xⁱᵗ[:,j]), w₀, w, v) for j in 1:sⁱ, i in 1:sᵘ]

# 事前計算あり
@time ŷᵘ = [y_hat(xᵘᵗ[:, i], 0, wᵘ, vᵘ) for i in 1:sᵘ]                              # 式(4)の第2項、第4項
@time qᵘ = precompute_q(xᵘ, vᵘ, sᵘ, k)                                            # 式(9)
@time ŷⁱ = [y_hat(xⁱᵗ[:, i], 0, wⁱ, vⁱ) for i in 1:sⁱ]                                  # 式(4)の第3項、第5項
@time qⁱ = precompute_q(xⁱ, vⁱ, sⁱ, k)                                               # 式(10)
@time e2 = [w₀ + ŷᵘ[i] + ŷⁱ[j] + qᵘ[:,i] ⋅ qⁱ[:,j] for j in 1:sⁱ, i in 1:sᵘ]  # 式(4)、式(8)

さて、実装のコツをふまえてコードを見ていきましょう。ベクトルや行列のサイズが記号定義と同じになっていますよね。また、簡単な内積演算は演算子を使っていますが、それ以外はfor文になっているので、疎行列の部分を除くと、ループの範囲が前節までの説明と一致していることも合わせて確認してみてください。

コンテキストデータ${\bf X}$を疎行列としているところと、行列にアクセスするときは列単位になっているところが実装上の注意点です。Juliaは基本的に列単位のアクセスが高速になるよう設計されているので、行単位でアクセスが必要な行列は転置して列単位でアクセスできるようにしています。${\bf Q}^u$や${\bf Q}^i$の定義で、因子を行、ユーザーやアイテムを列にしたのも列単位でアクセスしやすくするためです。

事前計算を行うことで少しコード量が増えていますが、この程度であれば実装工数の増加を懸念するほどではないでしょう。式($\ref{eq:fm3}$)の第2項と第4項は式($\ref{eq:fm}$)とほぼ同じ形なのでy_hat()関数を流用できます。式($\ref{eq:fm3}$)の第3項と第5項についても同様です。${\bf Q}^u$や${\bf Q}^i$の事前計算でコードが10行ほど増えますが、$\bf{Q}$同士の計算は内積計算なので簡単です。

最後に

今回は、FMの評価推定値の効率的な計算方法について紹介しました。技術的には大した方法ではないものの、少しの工夫で計算量を削減できるので実用性は高いです。このような実装の仕方は実務で数値計算やアルゴリズム実装などをされている方にとっては当たり前の方法ですが、短い工数で機械学習を使おうとすると計算量まで頭が回らないことがあります。そんなときに、こういう方法があることを知っていれば、役に立つこともあるのではないでしょうか。

以前の記事とともにFMを使ってレコメンデーションをするときの計算量に着目した内容を扱いました。機械学習では精度に注目が集まりがちでディープラーニング関係を除くと計算量の話題は少ないように思います。弊社の場合だと、グロース施策が当たってユーザー数が急増したり、業務提携でアイテム数(求人数)が急増したりすることがあるので、どのぐらいのデータ量までなら現在のままで対応できるか、インスタンス増強するとどれぐらいコストが増えるかなどの見通しをもっておく必要があります。こういうときに計算量は役に立ちます。機械学習は「やってみないとわからない」と言われますが、理解を深めることで見通しを立てられる場面も少なくありません。弊社では機械学習を使う場合でも見通しを立てたり原因を考えたりすることを大事にしています。

レコメンデーションは売上に直結することが多いのでよく使われますが、使いこなすのは大変です。RMSEやAUCのような精度に注目が集まりがちですが、今回紹介したように計算量が問題になることもあります。ログデータを使ってオフライン精度検証を行うときにも注意が必要で、単純にデータ分割する交差検証を使うのは不適切ですし、そのままRMSEを評価するとバイアスが生じる問題もあります。また、RMSEやAUCのような精度が高ければよいレコメンデーションになっているわけでもないので精度指標についても配慮が必要です。実際に運用してみないと気づきにくい課題もいろいろあるのですよね。そのような課題や対処方法は地味な内容だったりするので公開されていることは少ないですが実務では重要だったりします。本ブログでは、学習データの作り方や精度検証のコツといった他ではあまり扱われていない内容も今後取り上げていこうと思います。