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のような精度が高ければよいレコメンデーションになっているわけでもないので精度指標についても配慮が必要です。実際に運用してみないと気づきにくい課題もいろいろあるのですよね。そのような課題や対処方法は地味な内容だったりするので公開されていることは少ないですが実務では重要だったりします。本ブログでは、学習データの作り方や精度検証のコツといった他ではあまり扱われていない内容も今後取り上げていこうと思います。

Alternating Least SquaresによるFactorization Machinesのパラメータ推定

こんにちは、リブセンスで統計や機械学習関係の仕事をしている北原です。今回はレコメンデーションにも使えるFactorization Machines(FM)の効率的な学習アルゴリズムの紹介です。実装にはJuliaを使います。

実務で必要な要件を満たす機械学習ライブラリがなくて、機械学習モデルをカスタマイズすることってありますよね。最近はTensorFlowのような機械学習フレームワークが充実してきたので、そういう場合にはこれらのフレームワークを利用することが多いかもしれません。しかし、アルゴリズムによっては、フルスクラッチで実装することで大幅に効率化できるものもあります。今回扱うFMのAlternating Least Squares(ALS) はその一例です。そこで使われている効率化方法は実装が簡単でギブスサンプリングなどでも使うことができる便利なものなのですが、あまり知られていないようです。そこで、今回はこの計算効率化方法を紹介しようと思います。

Factorization Machines

レコメンドアルゴリズムといえば、協調フィルタリングやコンテンツベースフィルタリングが有名ですが、本記事ではコンテキストを使えるFMを扱います。求人サービスだと職種や年収、勤務地などの情報が有用なので、これらの情報をレコメンデーションで利用できると新規ユーザーにも好みにあったレコメンドがしやすくなります。FMは有名なので、検索するといろいろと情報が見つかると思います。ここでは計算効率化に関する部分を中心に説明します。

FMは、Matrix Factorization(MF)やSVD++のような協調フィルタリングを一般化して表現できるようにしたモデルです。MFと同様に、ユーザーが未評価のアイテムについても評価推定値$\hat{y}$を計算することができるので、レコメンデーションに利用されています。MFとFMは学習データのフォーマットが異なっていて、MFだと評価点数をユーザー数×アイテム数の疎行列で表したものを学習データとしますが、FMだと1行1評価としてユーザーやアイテムに対応する要素を1としたものを学習データとしているところに注意が必要です。詳しくは論文解説記事1解説記事2などを見てください。

FMは、評価データ数を$s$、ユーザーIDなども含むコンテキスト数(特徴量数)を$n$、因子数を$k$とし、評価値$\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}$を使って、以下のようなモデルで表現されます。

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

モデル上のポイントは、交互作用項の係数$\hat{w}_{i, j}$がモデルパラメータ${\bf V}$の異なる行の内積になっているところにあります。これにより特徴量の組み合わせが学習データに存在しなくとも$\hat{w}_{i, j}$を計算できるため、レコメンデーションに利用しやすくなっているというのが1つ目の利点です。通常の回帰モデルでは、交互作用項の係数が特徴量の組み合わせの数$n(n+1)/2$だけ必要です。そのため、特徴量数$n$が多いとモデルパラメータ数が膨大になるという問題があります。しかし、FMでは${\bf V}$のサイズ分しか必要としません。通常は$k \ll n$となるよう因子数$k$を決めるので、同様の問題は起きにくくなるというのが2つ目の利点です。

一見すると交互作用項が複雑なのですが、各モデルパラメータについては線形になっているところが計算上のポイントです。二乗誤差を考えたとき、各モデルパラメータについては単純な二次関数になるため、簡単に最適値を計算することができます。

ここで後の説明をやりやすくするために、式($\ref{eq:fm}$)の別の表記を導入します。各モデルパラメータについて線形なので、着目するモデルパラメータを$\theta$としたとき式($\ref{eq:fm}$)は一般的に

$$ \begin{equation} \hat{y}(\mathbf{x} | \theta)=g_{(\theta)}(\mathbf{x})+\theta h_{(\theta)}(\mathbf{x}) \label{eq:lin} \end{equation} $$

と書くことができます。$g_{(\theta)}(\mathbf{x})$は$\theta$を含まない部分、$h_{(\theta)}(\mathbf{x})$は$\theta$にかかっている部分です。$\theta$以外のモデルパラメータはすでに与えられていて定数になっていると考えるとわかりやすいかもしれません。少し具体例を見てみましょう。例えば、特徴量数$n=3$、因子数$k=2$のとき、式($\ref{eq:fm}$)は

$$ \begin{eqnarray} \hat{y}(\mathbf{x}) &=& w_0 + w_1 x_1 + w_2 x_2 + w_3 x_3 \nonumber \\ &+& (v_{1, 1} v_{2, 1} + v_{1, 2} v_{2, 2}) x_1 x_2 + (v_{1, 1} v_{3, 1} + v_{1, 2} v_{3, 2}) x_1 x_3 + (v_{2, 1} v_{3, 1} + v_{2, 2} v_{3, 2}) x_2 x_3 \end{eqnarray} $$

と書けます。$\theta=w_2$とすると

$$ \begin{eqnarray} g_{(w_2)}(\mathbf{x}) &=& w_0 + w_1 x_1 + w_3 x_3 \nonumber \\ &+& (v_{1, 1} v_{2, 1} + v_{1, 2} v_{2, 2}) x_1 x_2 + (v_{1, 1} v_{3, 1} + v_{1, 2} v_{3, 2}) x_1 x_3 + (v_{2, 1} v_{3, 1} + v_{2, 2} v_{3, 2}) x_2 x_3 \\ h_{(w_2)}(\mathbf{x}) &=& x_2 \end{eqnarray} $$

となります。同様に、$\theta=v_{3, 2}$とすると

$$ \begin{eqnarray} g_{(v_{3,2})}(\mathbf{x}) &=& w_0 + w_1 x_1 + w_2 x_2 + w_3 x_3 \nonumber \\ &+& (v_{1, 1} v_{2, 1} + v_{1, 2} v_{2, 2}) x_1 x_2 + v_{1, 1} v_{3, 1} x_1 x_3 + v_{2, 1} v_{3, 1} x_2 x_3 \\ h_{(v_{3,2})}(\mathbf{x}) &=& x_3 (v_{1, 2} x_1 + v_{2, 2} x_2) \end{eqnarray} $$

となります。式($\ref{eq:lin}$)の表記を使うことで、$w$や$v$の違いにかかわらず統一的にモデルパラメータを扱えたり、微分を

$$ \begin{equation} \frac{\partial \hat{y}(\mathbf{x})}{\partial \theta} = h_{(\theta)}(\mathbf{x}) \label{eq:par_y} \end{equation} $$

のように簡潔に表記できます。しかし、これだけだと大変な計算が$g_{(\theta)}(\mathbf{x})$や$h_{(\theta)}(\mathbf{x})$に置き換えられて数式が簡単になるぐらいで、実用上のメリットがないように思うかもしれません。特徴量数や因子数が増えれば$g_{(\theta)}(\mathbf{x})$や$h_{(\theta)}(\mathbf{x})$の計算量が増加することは容易に想像がつきます。ところが、$g_{(\theta)}(\mathbf{x})$や$h_{(\theta)}(\mathbf{x})$を定義どおりに計算しなくてもモデルパラメータを推定できるうまい方法があり、式($\ref{eq:lin}$)の形にすることでその方法を導出しやすくなっているのです。このうまい方法については後で説明します。

Alternating Least Squaresによるパラメータ推定

モデルパラメータの推定方法には様々なものがありますが、今回はAlternating Least Squres(ALS、交互最小二乗法)を使います。MFの派生モデルのパラメータ推定法としては定番の一つです。Coordinate Descent(座標降下法)と考えることもできます。最近だとStochastic Gradient Descent(SGD、確率的勾配降下法)のほうが有名ですし、FMのモデルパラメータ推定にも利用されていますが、SGDだと今回紹介する効率的計算は使えません。

ALSでは、更新しないモデルパラメータを固定したまま、誤差を最小化するように各パラメータを一つずつ更新する処理を、最適値に収束するまで繰り返します 。一方で、SGDやGradient Descent(GD、勾配降下法)のような勾配法では、勾配を少しずつ下っていくことで最適値をもとめます。WikipediaのCoordinate descentGradient descentに掲載されている図を見ると、直感的なイメージがつかめると思います。

SGDやGDのような単純な勾配法と比較したときのALSの利点は、学習率のような調整パラメータがないところです。近年では学習率を自動調整する方法がいろいろと提案されているので、あまり強いメリットを感じないかもしれません。それでも、最初から学習率を考えなくともよいのは利点の一つでしょう。

では、FMのモデルパラメータ更新式を導出しましょう。モデルパラメータ集合を$\Theta$、モデルパラメータ$\theta \in \Theta$の正則化係数を$\lambda_{(\theta)}$として正則化項を入れた二乗誤差

$$ \begin{equation} L = \sum^s_{l=1}{(\hat{y}({\bf x}_l) - y_l)^2} + \sum_{\theta \in \Theta}{\lambda_{(\theta)} \theta^2} \label{eq:err} \end{equation} $$

を最小化することでパラメータを推定します。モデルパラメータは複数ありますが、ALSでは1パラメータのみを変数として考え、残りを定数のように扱うところがポイントです。1つのモデルパラメータ$\theta$のみを変数としたとき、二乗誤差$L$を最小にする$\theta$の値は、$\theta$で微分したものを0とおいて整理すればよいですね。式($\ref{eq:lin}$)と式($\ref{eq:par_y}$)を使って

$$ \begin{eqnarray} \frac{\partial L}{\partial \theta} &=& \sum^s_{l=1}{2 (\hat{y}({\bf x}_l) - y_l) \frac{\partial \hat{y}({\bf x}_l)}{\partial \theta}} +2 \lambda_{(\theta)} \theta \nonumber \\ &=& 2 \sum^s_{l=1}{(\hat{y}({\bf x}_l) - y_l) h_{(\theta)}(\mathbf{x}_l)} +2 \lambda_{(\theta)} \theta \nonumber \\ &=& 2 \sum^s_{l=1}{(g_{(\theta)}({\bf x}_l) - y_l) h_{(\theta)}(\mathbf{x}_l)} + 2 \theta \sum^s_{l=1}{h^2_{(\theta)}(\mathbf{x}_l)} +2 \lambda_{(\theta)} \theta =0 \\ \end{eqnarray} $$

と変形し、$\theta$について整理すると

$$ \begin{equation} \theta=-\frac{\sum^s_{l=1}{(g_{(\theta)}({\bf x}_l) - y_l) h_{(\theta)}(\mathbf{x}_l)}}{\sum^s_{l=1}{h^2_{(\theta)}(\mathbf{x}_l)} + \lambda_{(\theta)}} \label{eq:orig_theta} \end{equation} $$

と、二乗誤差$L$を最小にする$\theta$を計算できます。式($\ref{eq:lin}$)を使っているのでシンプルな形にまとまっていますね。式($\ref{eq:fm}$)や式($\ref{eq:w}$)をそのまま愚直に微分してみると、違いが実感できると思います。

次に、1パラメータ1回の更新に必要な計算量について考えてみましょう。実際には全パラメータについて収束するまで値更新を繰り返すので、さらにパラメータ数と繰り返し数を乗じたものが全体の計算量になることに注意してください。

基本的には式($\ref{eq:orig_theta}$)の通りに計算すれば、$\theta$を算出できます。時間がかかるのは$\theta=v_{i, j}$のときなので、工夫なしに計算すると、$n$個のパラメータの組み合わせが$g_{(\theta)}(\mathbf{x}_{l})$の計算にあらわれるので計算量は$O(s \, k \, n^2 )$となります。レコメンデーションでは、コンテキスト数はユーザー数やアイテム数より少ないことが多いので、$n$はユーザー数とアイテム数の和と同じオーダーになることが多いです。問題設定にもよりますが弊社の場合、$n$のオーダーは数万から数十万程度になるので、$O(s \, k \, n^2)$だと現実的な時間で計算を終わらせるのは難しいことが多いです。

計算を効率化する方法はいくつかあって、FMが最初に提案された論文のLemma 3.1を利用すると、式($\ref{eq:fm}$)の交互作用項は

$$ \begin{equation} \sum^n_{i=1}\sum^n_{j=i+1}\sum^k_{f=1}{v_{i,f} v_{j,f}} x_i x_j = \frac{1}{2} \sum_{f=1}^{k}\left(\left(\sum_{i=1}^{n} v_{i, f} x_{i}\right)^{2}-\sum_{i=1}^{n} v_{i, f}^{2} x_{i}^{2}\right) \label{eq:interaction_term} \end{equation} $$

と変形できるので、計算量は$O(s \, k \, n)$になります。さらに、データがスパースな場合、1サンプルあたりの非ゼロデータの個数の平均を$\overline{m}_R$とすると、$O(s \, k \, \overline{m}_R)$になります。レコメンデーションなどで使う場合は、$\overline{m}_R \ll n$となることがほとんどです。例えば、ユーザー数とアイテム数が数万のオーダーとして$n \sim O( 10^4 )$、1サンプルあたりの非ゼロデータ数が百程度として$\overline{m}_R \sim O(10^2)$となっている場合、何も工夫をしなかったときと比較すると$O(10^6)$、つまり百万倍速くなります。これぐらいになると、現実的な時間内で計算が終わります。しかし、この計算を全パラメータについて行うと、交互作用項の計算に必要なパラメータ$\bf{V}$のサイズが$n \times k$ なので、$O(s \, n \, k^2 \, \overline{m}_R)$となり、因子数$k$を大きくすると計算時間が$O(k^2)$で急激に増えてしまいます。表現能力を高めるためにオーバーフィットしない範囲で因子数を大きくすることが多いですが、これだと気軽に大きくはできません。

事前計算による高速学習

さて、本題の効率的な計算方法について説明します。

計算を効率化するポイントは、共通する計算は一度で済ませ、差分のみ計算するところにあります。着目すべきは、全データサンプルについての和を計算している部分に含まれる、更新されるパラメータ$\theta$について線形な部分です。例えば、この線形な部分を$f(\mathbf{x}, \theta)=g(\mathbf{x}) + \theta h(\mathbf{x})$と書き、更新後のパラメータを$\theta^{\ast}$とすると$f(\mathbf{x}, \theta^{\ast})=g(\mathbf{x}) + \theta^{\ast} h(\mathbf{x}) = f(\mathbf{x}, \theta) + (\theta ^{\ast} - \theta) h(\mathbf{x})$のように更新後の$f(\mathbf{x}, \theta^{\ast})$を計算できます。$f(\mathbf{x}, \theta)$は更新前のものなので、前のステップで記憶しておけば計算量は$O(1)$です。そのため、$g(\mathbf{x})$の計算量がどれだけ多くとも、$h(\mathbf{x})$の計算量のみに依存することになります。前のステップで計算したものを使うので、ここではこの処理を事前計算と呼びます。

これをふまえ、効率化の方針としては、式($\ref{eq:orig_theta}$)からパラメータ$\theta$について線形な部分を探し、前のステップで記憶しておくべき部分と更新式の整理をすることで計算量の削減を図ります。式($\ref{eq:orig_theta}$)の分子は$g_{(\theta)}({\bf x}_l) - y_l$と$h_{(\theta)}(\mathbf{x}_l)$の積になっていて一度に扱うと複雑なので、各項について効率化を図っていきましょう。$g_{(\theta)}({\bf x}_l) - y_l$と$h_{(\theta)}(\mathbf{x}_l)$の計算量はそれぞれ$O(k \, \overline{m}_R)$、$O(\overline{m}_R)$なので、計算量の大きい$g_{(\theta)}({\bf x}_l) - y_l$から先に考えます。

残差$e(\mathbf{x}_l, y_l | \Theta) = \hat{y}(\mathbf{x}_l | \Theta) - y_l$を使うと

$$ \begin{eqnarray} g_{(\theta)}({\bf x}_l) - y_l &=& \hat{y}(\mathbf{x}_l | \Theta) - y_l - \theta h_{(\theta)}(\mathbf{x}_l) \nonumber \\ &=& e(\mathbf{x}_l, y_l | \Theta) - \theta h_{(\theta)}(\mathbf{x}_l) \end{eqnarray} $$

と書けます。$e(\mathbf{x}_l, y_l | \Theta)$は$\theta$について線形なので、更新後のパラメータを$\theta^{\ast}$を含むパラメータ集合を$\Theta^{\ast}$とすると

$$ \begin{equation} e(\mathbf{x}_l, y_l | \Theta^{\ast}) = e(\mathbf{x}_l, y_l | \Theta) + (\theta^{\ast} - \theta)h_{(\theta)}(\mathbf{x}_l) \label{eq:e} \end{equation} $$

と更新できます。更新前の$e(\mathbf{x}_l, y_l | \Theta)$を記憶しておけば、計算量は$O(\overline{m}_R)$になるので、$k$倍速くなりますね。

次に、$h_{(\theta)}(\mathbf{x}_l)$についても考えてみましょう。$\theta=w_i$の場合は、$h_{(w_i)}(\mathbf{x}_l) = x_{l, i}$なので常に$O(1)$となるので簡単です。$\theta=v_{i, f}$の場合、式($\ref{eq:fm}$)、式($\ref{eq:w}$)、式($\ref{eq:par_y}$)より

$$ \begin{eqnarray} h_{(v_{i, f})}(\mathbf{x}_l) &=& \frac{\partial \hat{y}(\mathbf{x}_l)}{\partial v_{i, f}} \nonumber \\ &=& x_{l, i} \sum_{j=1}^{n} v_{j, f} x_{l, j} - x^2_{l, i} v_{i, f} \label{eq:orig_h} \end{eqnarray} $$

と計算できます。交互作用項の微分は式($\ref{eq:interaction_term}$)の形を使ったほうが計算ミスしにくいかもしれません。ここで$q_{l, f}(\Theta) = \sum_{j=1}^{n} v_{j, f} x_{l, j}$に着目し$\theta=v_{j, f}$と考えると、$q_{l, f}(\Theta) $は$\theta$について線形なので

$$ \begin{equation} q_{l, f}(\Theta^{\ast}) = q_{l, f}(\Theta) + (v^{\ast}_{i, f} - v_{i, f}) x_{l, i} \label{eq:q} \end{equation} $$

と更新でき、

$$ \begin{equation} h_{(v_{i,f})}(\mathbf{x}_l) = x_{l, i} q_{l, f}(\Theta) - x^2_{l, i} v_{i, f} \label{eq:h} \end{equation} $$

と整理できます。つまり、$q_{l, f}(\Theta)$を記憶しておけば、$h_{(\theta)}(\mathbf{x}_l)$も$O(1)$で計算できるということです。そのため、式($\ref{eq:orig_theta}$)の分母にある$h^2_{(\theta)}(\mathbf{x}_l)$も$O(1)$で計算できるので、結局、式($\ref{eq:orig_theta}$)は$O(s)$で計算できることがわかります。

さらに、$x_{l, i}=0$だと、式($\ref{eq:q}$)、($\ref{eq:h}$)より$h_{\theta}(\mathbf{x}_l)=0$なので、計算の必要がありません。1パラメータに対応する非ゼロデータの個数の平均を$\overline{m}_C$とすると、式($\ref{eq:orig_theta}$)は$O(\overline{m}_C)$となります。レコメンデーションの場合、ユーザーIDやアイテムIDに対応するカラムはほとんどゼロで$\overline{m}_C \ll s$となることがほとんどなので、さらに計算量は少なくなります。

では、事前計算によってどのぐらい計算量が削減されているか考えてみましょう。事前計算を使わないときの計算量は$O(s \, k \, \overline{m}_R)$でしたから、事前計算によって$s \, k \, \overline{m}_R / \overline{m}_C$倍速くなります。通常は評価サンプルデータ数$s$のほうがコンテキスト数$n$より多いので、$\overline{m}_C < \overline{m}_R$となります。例えば、$s \sim O(10^5)$、$k \sim O(10^2)$、$\overline{m}_R \sim O(10^2)$、$\overline{m}_C \sim O(10)$とすると、一億倍速くなるということになります。もちろん、計算のオーバーヘッドなどがあるため実際にはこれほど速くなることは稀ですが、それでも計算量削減にかなり有効であることはわかると思います。

最終的に、パラメータ数のオーダーが$O(n \, k)$であることを思い出すと、全パラメータを一回ずつ更新するのに必要な計算量は$O(n \, k \, \overline{m}_C )$となります。これぐらい速いと、弊社の事業規模では問題になることがほぼありません。

実装

では実装してみましょう。ここではJuliaを使います。少し余談ですが、最近フルスクラッチでアルゴリズムを書く場合は簡潔な表記が可能で計算も高速なJuliaを使うことが多いです。別記事でも紹介しているように弊社には機械学習基盤があるので、使いたい言語やライブラリで実装したアルゴリズムを実サービスに導入しやすくなっています。ここで紹介するコードも社内で使われているものをベースにしています。

実装しやすいようにモデルパラメータ$w_0$、$\bf{w}$、$\bf{V}$ごとに式($\ref{eq:orig_theta}$)を整理すると

$$ \begin{eqnarray} w_0 &=& -\frac{\sum^s_{l=1}{e({\bf x}_l, y_l | \Theta) - s w_0}}{s + \lambda_{(w_0)}} \\ w_j &=& -\frac{\sum^s_{l=1}{(e({\bf x}_l, y_l | \Theta) - w_j x_{l, j}) x_{l, j}}}{\sum^s_{l=1} x^2_{l, j} + \lambda_{(w_j)}} \\ v_{j, f} &=& -\frac{\sum^s_{l=1}{(e({\bf x}_l, y_l | \Theta) - v_{j, f} h_{(v_{j, f})}({\bf x}_l)) h_{(v_{j, f})}({\bf x}_l)}}{\sum^s_{l=1}h_{(v_{j, f})}({\bf x}_l) + \lambda_{(v_{j, f})}} \end{eqnarray} $$

となります。 これをコードに落とすと、モデルパラメータ推定のコードは以下のようになります。

using Distributions
using SparseArrays
using LinearAlgebra

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²  # 式(13)
    end
    return w₀ + x ⋅ w[i] + 0.5 * ΣΣwᵢⱼxᵢxⱼ
end

function precompute(y, x, w₀, w, v)
    s = size(x)[1]
    k = size(v)[2]
    e = zeros(s)     # 残差
    q = zeros(s, k)  # 交互作用項のカラム非依存部分
    xᵗ = sparse(x')  # CSC行列なので行単位でアクセスする場合は転置して列単位でアクセス

    for l in 1:s
        e[l] = y_hat(xᵗ[:, l], w₀, w, v) - y[l]
        for f in 1:k
            (l_, x_) = findnz(xᵗ[:, l])
            q[l, f] = x_ ⋅ v[l_, f]
        end
    end
    return [e, q]
end

function train(y, x, k, λ, n_itr)
    s = size(x)[1]  # データ数
    n = size(x)[2]  # コンテキスト数(特徴量数)

    # Initialize the model parameters
    w₀ = 0.0
    w = zeros(n)
    v = rand(Normal(0.0, 0.001), n, k)

    # 事前計算
    e, q = precompute(y, x, w₀, w, v)
    h = zeros(Float64, s)

    # main optimization
    for itr in 1:n_itr
        # global bias
        w₀ᵃ = -(sum(e) - s * w₀) / s
        e .+= (w₀ᵃ - w₀)  # 式(15)
        w₀ = w₀ᵃ

        # 1-way interactions
        for j in 1:n
            (i_, x_) = findnz(x[:, j])  # 非ゼロデータのインデックスと値

            wᵃ = -sum((e[i_] .- w[j] .* x_) .* x_) / (sum(x_.^2) + λ)  # 式(20)
            e[i_] .+= (wᵃ - w[j]) .* x_                                # 式(15)
            w[j] = wᵃ
        end

        # 2-way interactions
        for f in 1:k
            for j in 1:n
                (i_, x_) = findnz(x[:, j])  # 非ゼロデータのインデックスと値

                h[i_] .= x_ .* q[i_, f] .- x_.^2 .* v[j, f]
                vᵃ = -sum((e[i_] .- v[j, f] .* h[i_]) .* h[i_]) / (sum(h[i_].^2) + λ)  # 式(21)
                e[i_] .+= (vᵃ - v[j, f]) .* h[i_]                                      # 式(15)
                q[i_, f] .+= (vᵃ .- v[j, f]) .* x_                                     # 式(17)
                v[j, f] = vᵃ
            end
        end
    end
    return [w₀, w, v]
end

Juliaを使っていることもあり、計算式をほぼそのままコードに置き換えるだけで実装できていることがわかると思います。実装上のポイントは、疎行列を使っているところと、その疎行列にカラム単位でアクセスしているところです。Juliaの疎行列はCompressed Sparse Column(CSC)形式なので、カラム単位のアクセスが高速です。そのため、事前計算で行単位のアクセスが必要なところでは転置した行列を使っています。ALSでは元々カラム単位で処理するので、パラメータ更新時に使われるデータは転置していません。

なお、コードをシンプルにするため、今回のテーマから外れる事項については、簡単な方法を使っています。例えば、ALSでは、モデルパラメータの更新順序が推定値の収束速度に影響しますが、本記事の主題から外れるので、オリジナル論文と同じく添え字順に$w_i$を更新後、$v_{i, j}$を更新する方法になっています。また、計算停止条件も固定回数で停止させる簡単な方法を使っています。正則化係数$\lambda_{(\theta)} $もモデルパラメータによらず一定の$\lambda$にしています。このあたりの話題も機会があれば記事にするかもしれません。

呼び出しコードは以下のようになります。評価数や因子数を増やしても線形にしか計算時間が増えないことを確認してみてください。

# サンプルデータ作成
s = 100           # 評価数
n = 100           # コンテキスト数
sparseness = 0.5  # 非ゼロの比率
y = rand(1:5, s)  # 1~5点の評価がついている想定
x = convert(SparseMatrixCSC{Int64,Int64}, sprand(Bool, s, n, sparseness))
xᵗ = sparse(x')   # CSC行列なので列単位でアクセスできるよう転置

# 調整パラメータ
k = 10     # 因子数
λ = 0.1    # 正則化係数
n_itr = 5  # 何回更新するか
@time w₀, w, v = train(y, x, k, λ, n_itr)

# 誤差が小さくなっているか確認
e = [y_hat(xᵗ[:, l], w₀, w, v) - y[l] for l in 1:n]
rmse = e.^2 |> mean |> sqrt

# 評価数を10倍にすると、計算時間もほぼ10倍になる
n = 1000
k = 10
x = convert(SparseMatrixCSC{Int64,Int64}, sprand(Bool, s, n, sparseness))
@time w₀, w, v = train(y, x, k, λ, n_itr)

# 因子数を10倍にしても、計算時間は10倍にしかならない
n = 100
k = 100
x = convert(SparseMatrixCSC{Int64,Int64}, sprand(Bool, s, n, sparseness))
@time w₀, w, v = train(y, x, k, λ, n_itr)

最後に

今回は、効率的な学習が可能になるALSによるFMのパラメータ推定について紹介しました。少しの工夫で計算量を大幅に削減しているところが実用的です。FMはMF派生モデルの一般化になっているというところに目がいきがちなのですが、計算量の点でも優れていることがわかると思います。この計算方法は弊社の一部のサービスで利用されていますが、モデルパラメータ推定の計算時間が問題になったことは今のところありません。計算量が問題になっていると並列・分散処理で何とかしようと考えてしまいがちですが、力技だけだとコストがかかりすぎる場合はこのような方法も検討する価値があるかもしれません。

今回はFMのモデルパラメータ推定を扱いましたが、FMをレコメンデーションで使う場合は評価推定値計算の計算量も問題になりやすいです。なぜかというと、各ユーザーに対してレコメンド対象アイテムの評価推定値を計算するので、組み合わせが多くなってしまうためです。学習データで交差検証を行ったりコンペで使ったりしただけだと問題にならないので気づきにくいのですよね。この話題については、後日、別の記事で扱う予定です。

DigdagとEmbulkで行うDB同期の管理

データプラットフォームグループの松原です。 弊社各サービスのデータ分析基盤であるLivesense Analytics(以降LA)の開発、運用を行っています。
今回はLAで行っている分析のためにサービス側のデータ(テーブル)を、Redshiftへ同期を行う処理について紹介します。

概要

LAではデータウェアハウスとしてRedshiftを運用しており、社内から比較的自由に利用できる様にしています。
LAで取り扱っているデータはアクセスログが中心ですが、分析を行う利用者からはLA由来のデータ以外にも自分たちのサービスのデータを用いて分析を行いたい、という要望がよく出てきます。
サービスのデータには個人情報を含むものも少なくありませんが、分析基盤として社内にデータを解放するためにはそのような情報は削る必要があります。
そこで個人情報をマスキングしたサービス側データを利用できるよう、Redshiftに同期しています。

システム構成(概要)

大まかなシステム構成としては、下図のようになっています。

f:id:livesense-analytics:20190227151830p:plain:w650

多くのテーブルを同期しているので、カラム単位のマスキング処理や加工処理はEmbulkで行わず、同期元のデータベースであらかじめ行っておきます。

Digdagの選定理由

BigData-JAWS勉強会でAirflowのことを話してきましたにもあるように、すでにAirflowを運用していますが、ここではDigdagを選定しています。
主な選定理由としては以下があります。

  1. 構築や運用の容易さ
    Digdagはサーバーモードで動かしても、JavaとPostgreSQLの動く環境があれば運用を開始できます。
    また、Airflowと比べると、今のところ構築・アップグレードなどの運用は容易です。
  2. やりたいことがシンプル
    今回はsh(embluk)、redshift、httpオペレータで用途として足りそうで、それぞれのタスク間の依存関係も複雑では無いです。
    そのためタスク間の依存関係はYAMLで宣言的に書いたほうがメンテンスしやすいと考えました。
    また、DAG自体も複雑で無いため、どこまで実行できたか・どこで失敗したのかの可視性は求めていないのでダッシュボードも簡易なもので問題ありませんでした。

これらの理由により、Airflowで統一するより用途によりツールを使い分けた方が運用が容易なのでは、という結論になりました。

規模感

Redshiftへ同期しているデータに関する情報ですが以下のとおりです。

同期元のデータベース: 13個
同期を行っているテーブル数: 280テーブル
同期対象のデータベースの種類: 2種類(MySQL、PostgreSQL)
同期先: 2種類(Reshift、Redshift Spectrum)

これらのテーブルは、毎日業務開始前に前日のデータを同期しています。
同期はデータの増分を追加するのではなく、全量差し替えています。

Digファイルの依存関係

現時点では同期対象のテーブル数が280テーブル程度あるため、テーブル個別のデータ変換処理はなるべく行わないようにしています。
embulkやdigファイルを共通化しておくことで、同期対象テーブルの追加などの依頼があった場合にはテーブルとカラムの情報が書かれたファイル(下記の例だとtest1_tables.dig)のみを変更すれば良いようにしています。
また、テーブル個別のデータ変換処理はできるだけ行わないようにしていますが、変換が必要な場合はAirflow側で処理を行うようにしています。

2つのDBと同期を行う際のdigファイルの依存関係は、下図ようになります。

それぞれのファイルは以下のような役割で分けています。

ファイル名 概要
test1.dig Workflowを表すdig
test1_secret.dig 同期元への接続情報を保持しているdig
test1_tables.dig 同期対象のTable,Columnを定義しているdig
from-mysql-to-redshift.dig Embulkの実行・Redshiftへの処理を定義しているdig
in-mysql-out-s3.yml.liquid Embulkの設定ファイル
create.sql 一意になるテーブル名を作成
copy.sql EmbulkでS3においたファイルを、create.sqlで作成したテーブルにCOPY
swap.sql 既存のテーブルとSWAPし既存のテーブルを削除

スケジュール設定、通知・エラー処理などを省略すると、主だったdigファイルは以下のようになります。

--test1.dig
_export:
  mysql_database: MySQLの接続先DB名(ex. tes1)
  la_schema: redshift上のスキーマ名
  !include : digdag/environment/env.dig
  redshift:
    host: redshiftのホスト名
    database: redshiftのDB名
    user: la_importer
    ssl: true
    schema: ${la_schema}
    strict_transaction: false
    connect_timeout: 600s

+task:
  _export:
    !include : digdag/secret/test1_secret.dig

  for_each>:
    !include : digdag/media/test1_tables.dig
  _do:
    !include : digdag/flow/from-mysql-to-redshift.dig
--digdag/media/test1_tables.dig
mysql_table_info:
  - table: table1
    column: id, created_at, updated_at
  - table: table2
    column: id, created_at, updated_at
--digdag/flow/from-mysql-to-redshift.dig
_export:
  la_table: ${mysql_table_info.table}
  s3_file_path: mediadb/${la_schema}/${la_table}

+from-mysql-to-s3:
  sh>: embulk run -b /home/digdag/embulk_bundle /home/digdag/dag/embulk/in-mysql-out-s3.yml.liquid
  _export:
    mysql_table: ${mysql_table_info.table}
    mysql_select: ${mysql_table_info.column}

+mysql_redshift_create_temptable:
  redshift>: digdag/flow/queries/create.sql

+mysql_redshift_copy_2_temptable:
  redshift>: digdag/flow/queries/copy.sql

+mysql_redshift_swap_table:
  redshift>: digdag/flow/queries/swap.sql

運用上行っていること

以下のようなことを行って、運用コストを下げています。

  • Digdagで同期するテーブル情報の作成・Redshift用のCreate Table文(SORTKEY, DISTKEYは手動で設定)の作成のスクリプト化
  • タスクのログをS3に出力し、(shオペレーターから呼び出した)EmblukのログをLambdaでパース処理を行いエラー内容をSlackへ通知
  • Redashで同期済みデータのヘルスチェック
  • Errbotとmogを用いてChatOps(DAGの実行)
  • Redshift上に同期したが、1年以上利用されてないテーブルの割り出しと、今後も同期するのかの確認(不要テーブルを同期対象から外す)

これらを含めると、全体では以下のような構成になります。
f:id:livesense-analytics:20190220113717p:plain:w650

まとめ

もともとは Airflow を用いたデータフロー分散処理にあるバッチ処理と同様にRakeとwhenever(Cron)で運用されていたものをDigdagに置き換えました。
置き換えを行うことで、メンテナンスコストはだいぶ下げることが出来たのではないかと思います。
現在はデータソースがどこにあるのかと処理の複雑さによって、AirflowとDigdagのどちらを使うのかを使い分けています。

今後やっていきたいこと

今後もDigdagを用いた同期の改善の他にも、以下のようなこともやってきたいと考えています。

  • DigdagとAirflowの連携の強化(同期完了をトリガーとしてAirflowのDAG実行)
  • マスキング処理にProxySQLを使うことで、任意のタイミングでの同期の実行
  • CDC(Change Data Capture)を用いた変更内容の(ほぼ)リアルタイムの同期

UXデザインが総論賛成、各論疑問になる理由と、プロジェクト設計で意識したい3つの条件【中編】

 テクノロジカルマーケティング部 データマーケティンググループにてUXリサーチャーをしている佐々木と申します。普段は、UXデザイン(以下、UXDと略記)に関するプロジェクトを事業部横断で支援する業務についております。

 前回の私のブログでは、前編として「UXデザインが総論賛成、各論疑問になる理由」と「プロジェク設計で意識したい3つの条件の<その1>サービス・ドミナント・ロジックを中核とするサービスかどうか」について述べさせて頂きました。 analytics.livesense.co.jp

今回は、中編として<その2>潜在ニーズを探る必要性の有無、について述べさせて頂きます。

前編でも述べさせて頂きましたがUXDのプロジェクト設計で意識したい3つの条件は

  • <その1>サービス・ドミナント・ロジックを中核とするサービスかどうか
  • <その2>潜在ニーズを探る必要性の有無
  • <その3>改善・改良を目的とした既存事業・既存サービス

になります。今回は3条件の中の2つ目です。

4.1 潜在ニーズ有無による特徴

 プロジェクト設計で意識したい2つ目の条件として、「潜在ニーズを探る必要の有無」を挙げさせて頂きました。 ここでは、まず、潜在ニーズの有無による新規ソリューション案(事業案)の特徴について考えてみます。

 一般的に「顕在ニーズ」を対象にしたソリューション案(事業案)は、購買欲求が高く自明な顕在ニーズを対象にしているので早期に市場が立ち上がることが期待できます。しかしながら、デメリットとして、競合が多く、差別化を作り難く、価格競争になりがちな傾向があります。

 一方、「潜在ニーズ」を対象にしたソリューション案(事業案)は、近未来における社会やユーザーの変化に対応したサービスにより差別化を作り易く、先行利益を得やすいです。しかしながら、デメリットとして、サービス市場が立ち上がるまで安定した需要が見込めず事業成長の不確実性が高くなる傾向があります。

 この様な特徴を考慮して「潜在ニーズを探る必要の有無」により、そのソリューション案(事業案)の種を探す工数や、説得材料を揃えて検証する工数を大きく変えて設計した方が良いと考えてます。図4.1.1に示したデザインプロセスにおけるダブルダイヤモンドで例えると、前後のフェーズで共に設計を変える必要があるのです。

f:id:livesense-analytics:20180926172204p:plain
図4.1.1 デザインプロセスにおけるダブルダイヤモンド

4.2「問題の発見・収束」における設計の違い

 ダブルダイヤモンドの前半である「問題の発見・収束」のステップについて、新規事業や新規サービスを検討する場合について考えます。対象となる問題が「顕在ニーズ」「潜在ニーズ」のどちらなのか?と、その新規性が「機能」もしくは「体験価値」の場合で考えてみます。

4.2.1 「顕在ニーズ」の場合

 差別化要素となる新規性が機能に有り、その機能が解決するのが「顕在ニーズ」の場合、どの様な問題(顕在ニーズ)を対象にするのかが自明なのでダブルダイヤモンドの前半である「問題の発見・収束」フェーズに掛ける工数は少なくすることができます。

 潜在ニーズを探る必要が余り無いのは、検討しているサービスの事業領域において世の中のサービスに不満が多く存在し”なんとなく不満なサービス”が世の中に溢れている状況になっている時、とも言えます。

 例えば図4.2.1 の様に「事業やサービス領域において、既存技術や新規技術で解決できる顕在ニーズが一定規模以上で存在している」状況などを挙げることができます。 ※図4.2.1 の丸角四角は、個別の顕在ニーズを表し、その大きさはニーズを持つ人の規模を表現しています。青塗りの丸角四角は、すでにソリューションが世の中に存在しているニーズ、白塗りの丸角四角は、まだソリューションが世の中に存在していないニーズ、を表します。

 高度成長時代に技術イノベーションにより産業成長が続いていた時は、【Ⅰ】【Ⅱ】を繰り返していたとも言えます。ですので比較的に技術的な発展による新興サービスが該当することが多いです。また、既存事業等で解決したい問題が明らかになっており、その解決案の目処がたっている時も該当すると言えます。

 この様な状況下では、従来の「定量的なマーケットリサーチの判断基準」で判断することができます。ダブルダイヤモンドの左側でユーザーのニーズについて深いリサーチをかける必要は余りありません。定量的に把握可能な自明な顕在ニーズに対して、素早く製品・サービス開発を行い市場に投入した方が事業成果を得やすいことが多いからです。

f:id:livesense-analytics:20180531115416p:plain
図4.2.1 潜在ニーズを探る必要が余り無い状況

4.2.2 「潜在ニーズ」の場合

 差別化要素となる新規性が機能で有るものの、機能により解決される一定規模感の有る顕在ニーズが見当たらない場合や、機能に差別化する要素が乏しく差別化要素となる新規性を体験価値に求める場合は、潜在ニーズの探索が必要な場合が多いです。また、自明になっていない問題(潜在ニーズ)を対象にするので「問題の発見・収束」フェーズに掛ける工数を多く設計する必要が有ります。

 潜在ニーズを探る必要のあるのは、事業環境面では検討しているサービスの事業領域において世の中のサービスに大きな不満は無く、”なんとなく良いサービス”が世の中に溢れている状況になっている時、とも言えます。

 例えば図4.2.2の【III】【Ⅳ】の様に「事業やサービス領域において、既存技術で解決できるニーズはほぼソリューションが存在しており」、【III】「既存技術で解決できる未解決の顕在的なニーズは少数意見のみ」の状況や、【Ⅳ】は【III】に加えて「新規技術で解決できる未解決の顕在ニーズも少数意見のみ」の状況を挙げることができます。

 この様な状況は、成熟期の事業領域に多く、従来のマーケットリサーチの判断基準(一定数量のある未解決なニーズを対象にサービス開発を行う)だけでは、新規事業の芽を見つけ出すことは非常に困難です。この時、様々なアプローチ方法がありますが、その一つに潜在ニーズを探る手段があります。

 潜在ニーズを探す場合、ダブルダイヤモンドの左側で、"今まで、顕在ニーズとして少数の人しか必要としていないと思われていたけど、実は多くの人が求めていた潜在ニーズ"を探しだすことを目的とします。

 ここで述べる潜在ニーズとは、厳密な定義とは少しずれてしまうかもしれませんが、ウォンツとニーズにおける”ウォンツ”、ジョブ理論における”本質的に片付けたいジョブ”、UXDにおける”本質的に求めている未発見の体験価値”、等と重なる概念として考えております。

f:id:livesense-analytics:20180531115653p:plain
図4.2.2 潜在ニーズを探る必要のある状況

4.3 「解決案の発見・収束」における設計の違い

 次にダブルダイヤモンドの後半である「解決案の発見・収束」のステップにおいて考えてみます。「解決案の発見・収束」のステップでは、前半で対象とした「顕在ニーズ」「潜在ニーズ」のどちらの解決策を探索するのかにより必要な工数は異なります。

 「顕在ニーズ」を対象にした解決策は、顕在ニーズを対象にしているので想定価格や購入意向率を客観的に評価することができます。よって評価と検証がしやすいことから工数を少なくすることができます。

 しかしながら「潜在ニーズ」を対象にした解決策は、まだ無消費の市場であることが多いため客観的な評価はありえません。解決策はステップを経た検証が必要です。またその説得材料を揃えるにも不確実性が高く限られた材料になるので時間が掛かります。数ターンのアジャイルなプロトタイプの検証を通して、想定価格や購入意向率等を、その「世界観」や「体験価値を実感できるストーリー」のプロトタイプを用いて検証し徐々に確度を上げていく必要があります。 またその際は、不確実性が高い検証になりますので、有力な一つの「潜在ニーズを対象にした解決策」に絞るよりも、複数案を対象に検証す進め、その中から選別するプロセスを推奨します。

 更に最終的な検証はサービスのリリース後も継続されます。将来起こるであろう「社会の行動や価値感の変容」と共にサービスが世の中に浸透することになるからです。ですので、今回のプロジェクトで、どこまでの範囲を検証対象とするのか?を事前に確認して設計する必要が有ります。

 また複数の候補案が有り、「顕在ニーズ」「潜在ニーズ」が混在している場合は特に注意が必要です。同じ設計で同時並行的に進めようとすると「潜在ニーズ」を対象にした検証が時間不足で不十分になりがちで、中途半端な結果となり最終候補案として見送られてしまうことが多いからです。その必要性により十分な工数を確保して検証することを推奨します。

 この様に、事業としての成果を出すためには、一見同じ様なテーマのリサーチプロジェクトでも、「潜在ニーズ探索の必要性有無」により、大きく設計を変える必要があると考えてます。

5. 中編のまとめ

 以上、UXDのプロジェクトを選定する時に意識する3条件の2つ目である <その2>潜在ニーズを探る必要性の有無、について実践時に心掛けてきたことを中心にまとめてみました。

 次回は、後編として最後の条件である

<その3>改善・改良を目的とした既存事業・既存サービス

について述べさせて頂きます。

階層ベイズによる小標本データの比率の推定

こんにちは、リブセンスで統計や機械学習関係の仕事をしている北原です。今回は階層ベイズを使った小技の紹介です。推定にはStanを使います。

Webサービスに限らないかもしれませんが、CVRやCTRなど比率データを扱うことって多いですよね。弊社の求人サービスは成果報酬型であるため、各求人の採用率などを知りたいこともよくあります。しかし、求人別だとバイト求人や転職求人の応募数はそれほど多くないので、採用数を応募数で単純に割っただけでは極端な採用率になりがちです。今回は、このような分母の値が小さい比率のデータを、階層ベイズを使って計算する方法を紹介します。


応募数が少ないときの採用率計算の問題

まず、応募数が少ない求人の採用率計算が必要な理由と、このようなサンプルサイズが小さいデータの比率計算の問題について説明します。
その問題をふまえて、今回どのような推定を行いたいかを説明します。

弊社の成果報酬型求人サービスでは、応募だけでなく採用してもらうことで売上が発生します。そのため、応募率だけでなく採用率も重要になってきます。例えば、あまりに採用率が低い求人があれば、その求人にマッチした求職者に応募してもらうような施策を考えたり、採用方針を企業に問い合わせたりする必要があります。また、採用率が高い求人があれば、その理由を調べて他の求人の採用率を上げる施策に活かすこともできます。

採用率は採用数を応募数で割ることで計算しますが、バイト求人にしても転職求人にしても、求人ごとに集計すると安定した比率を得るだけの十分な応募数がある求人は少ないです。集計期間を長くとれば計算に使える応募数が多くなることもありますが、求人は入れ替わることも多いし、集計期間を長くすればするほど古い情報を反映した採用率になってしまうので好ましくありません。また、応募数が少なかった求人についても状況に応じて適切な対策を打ちたいので、応募数が少ない求人を無視するということもしたくありません。

当然ながら、採用率を計算するときの分母にあたる応募数が少ないと、算出された採用率の信頼性は低くなります。例えば、応募数200件で採用数100件の求人であれば採用率50%と言われてもあまり違和感がないですが、応募数2件で採用数1件だと採用率50%と言われても疑わしいですよね。特に応募数が少ないときに、極端な採用率が算出されてしまいやすいという問題があります。

それでは、全求人の平均採用率が20%であることがわかっていたらどうでしょう。多くの人は、各求人の応募数や採用数だけでなく全体的な傾向も考慮して、応募数が増えたときの採用率を推測するのではないでしょうか。例えば、応募数2件で採用数1件の求人の採用率は、応募数が増えていけば50%より低くなりそうですが、全体採用率20%よりは高そうなので、だいたい30〜40%ぐらいかなと考えませんか。一方、応募数200件で採用数100件の求人場合だと、分母の値も比較的大きいので採用率は50%より少し低めの45〜49%ぐらいかなと考えませんか。

このように全体的な傾向も考慮して採用率を計算できると便利です。今回はベイズ推定を使って、このような計算を行います。

主観的な事前分布を設定したベイズ推定

それでは、全求人の採用率傾向を考慮した各求人の採用率推定の考え方について説明します。

ベイズ推定では主観的な情報を事前分布として設定することができます。そのため、全求人の採用率の傾向を事前分布として設定すれば、全体傾向を考慮して各求人の採用率を推定することができそうです。

ベイズ推定では、データが少ないと事前分布の影響を強く受けますが、データが多くなってくると徐々に事前分布の影響が薄れてきます。そのため、応募数が少ないと、推定採用率はデータから単純に計算された採用率よりも全求人の採用率に近い値になります。一方で、応募数が多くなると、全求人の採用率よりデータから計算された採用率に近い値になります。その結果、前節で例として説明したような、人が直感的に判断したときと同じような推定採用率が得られることが期待できます。

事前分布の検討

今回の肝となる事前分布をどうするかについて説明します。

どのような事前分布を設定するのがよいかについては様々な考え方があるのですが、ここでは全求人の採用率傾向を表すことができ、計算が収束する分布を選ぶことにします。全求人の採用率傾向についてはデータを調べて検討することにしましょう。

ここからはサンプルデータを使って具体的に見ていきます。

サンプルデータ(クリックすると展開されます)


このサンプルデータは、ある条件で抽出した求人の応募数と採用数を模擬したデータになっています。1行が1求人に対応し、N、R、G列は推定に使うデータ、idとrateは参考のために付与したデータになっています。Nは応募数、Rは採用数、Gは地域や職種など何らかのグループを表すものと考えてください。rateは採用率で、idは説明の際に参照しやすいようにつけました。全体の採用数を応募数で割り、全体採用率を計算するとおよそ0.1766になります。全391求人のうち、応募が10件以上ある求人でも53しかなく、ほとんどの求人では採用率を計算しにくいことがわかると思います。

全体傾向を事前分布にするためには、採用率がどのような確率分布になるかを考える必要があります。まずは各求人の採用率をヒストグラムにして分布を見てみましょう。
f:id:livesense-analytics:20181004181021p:plain:w360
採用率0%が飛び抜けて多く、その他33%、50%、100%が多いことがわかりますが、これだけだと応募数が増えたときにどのような確率分布になるかわからないですね。知りたいのは応募数が多くなったときの採用率の分布なので、応募数が一定数以上の求人の採用率分布を見てみましょう。
f:id:livesense-analytics:20181004181227p:plain:w360f:id:livesense-analytics:20181004181239p:plain:w360f:id:livesense-analytics:20181004181243p:plain:w360
応募数が多い求人だけを見てみると、右裾の長い分布になっているように思えます。加えて、採用率0%や100%の求人もなくなっていることから、応募数が十分多いと採用率0%や100%の求人は0件になるのではと予想されます。そこで、このような特徴をもつ確率分布を探すとベータ分布が見つかります。これ以降はベータ分布を事前分布として使うことにしましょう。なお、採用率推定の事前分布としてベータ分布が常によいわけではありません。例えば、応募数が多いのに採用率100%の求人が多い場合はガンマ分布など別の分布を使うこともあります。

事前分布を活用した比率の推定

事前分布に使う確率分布も決めたので、実際に採用率を推定します。ここでは、前節で検討した事前分布の有無による違いを見るため、事前分布を使わないときと、事前分布に固定パラメータのベータ分布を使ったときの2種類の結果を調べます。

まず、採用率は二項分布にしたがうことだけを仮定したシンプルなモデルで推定してみましょう。Stanでは明示的に事前分布を指定しないと一様分布を指定したときと同じになります。モデルとコードは以下のようになります。
f:id:livesense-analytics:20181004182305p:plain:w480

Rコード(クリックすると展開されます)

library(rstan)
library(tidyverse)

dat <- read_csv("sample.csv")

stan_dat <-
  list(M = nrow(dat),
       N = dat$N,
       R = dat$R)

fit <- stan(file = 'binom.stan', 
            data = stan_dat,
            seed = 123)


Stanコード(binom.stan)(クリックすると展開されます)

data {
  int M;
  int<lower=0> N[M];  // 応募数(分母)
  int<lower=0> R[M];  // 採用数(分子)
}
 
parameters {
  vector<lower=0, upper=1>[M] theta;  // 採用率
}
 
model {
  R ~ binomial(N, theta);
}


推定結果のサマリは以下のようになります。
推定結果のサマリ(クリックすると展開されます)

単純に計算したときのように0%や100%のような推定採用率にはならないことがわかります。一方で、採用率の全体傾向がまったく考慮されていないため、応募1件採用0件の求人の採用率期待値がおよそ33%になるなど、採用率としては好ましくない結果も見られます。

では、採用率の事前分布としてベータ分布を設定したモデルで推定してみましょう。ここではベータ分布のパラメータに\(\alpha=3\)と\(\beta=13\)を使います(理由は後述します)。固定パラメータを使うのでプレート図は同じで、モデルとコードは以下のように事前分布の部分が追加されます。

f:id:livesense-analytics:20181004182611p:plain:w480

Rコード(クリックすると展開されます)

library(rstan)
library(tidyverse)

dat <- read_csv("sample.csv")

stan_dat <-
  list(M = nrow(dat),
       N = dat$N,
       R = dat$R)

fit <- stan(file = 'binom_beta.stan', 
            data = stan_dat,
            seed = 123)


Stanコード(binom_beta.stan)(クリックすると展開されます)

data {
  int M;
  int<lower=0> R[M];  // 応募数(分母)
  int<lower=0> N[M];  // 採用数(分子)
}
 
parameters {
  vector<lower=0, upper=1>[M] theta;  // 採用率
}
 
model {
  R ~ binomial(N, theta);
  theta ~ beta(3, 13);     // 事前分布
}


推定結果のサマリは以下のようになります。
推定結果のサマリ(クリックすると展開されます)


まず、採用がないケースについて見てみると、応募1件採用0件の求人は採用率期待値がおよそ0.176、応募2件採用0件の求人は採用率期待値がおよそ0.166となっています。応募数が1〜2件程度と少なく不確実性が高いので、全体採用率0.1766よりやや低い値として算出されています。また、応募1件採用0件より応募2件採用0件のほうが採用率はより低い可能性が高いため、採用率も低めになっていることが確認できます。さらに、id=234の求人のように応募10件採用0件になると、採用率が低い可能性がさらに高まるので、採用率期待値はおよそ0.116と推定されています。逆に採用が多いケースについて見てみると、応募1件採用1件の求人の採用率期待値はおよそ0.235、応募3件採用3件の求人ではおよそ0.314と推定されています。データが少なく不確実性が高いことを考えると、これらの推定結果は直感的によさそうに思えます。

応募が多いほど推定もより正確になるので採用率の分布の幅はより狭くなります。例として、単純に計算した時の採用率が等しいid=302とid=308の二つの求人について見てみましょう。id=302は応募4件採用1件、id=308は応募72件採用18件の求人です。推定された各々の求人の採用率の分布は以下のようになります。

f:id:livesense-analytics:20181004183102p:plain:w360

id=302は応募数が少ないので事前分布の影響を強く受けており、さらに不確実性も高いので分布の幅がid=308と比較して広いことがわかります。一方で、id=308は応募が比較的多いため事前分布の影響が弱く、さらに不確実性も低いので分布の幅が狭いことがわかります。

目論見どおりベータ分布を事前分布として設定することで、直感的によさそうな採用率が推定されているように見えます。しかし実は、推定結果がベータ分布のパラメータ\(\alpha\)、\(\beta\)の値に強く依存するという問題があります。そのため、ベータ分布のパラメータ\(\alpha\)、\(\beta\)を別の固定値に設定すると、全く異なる推定採用率になってしまいます。採用率のヒストグラムを見ても、応募の多い求人が十分にないため適切なパラメータ値まではよくわかりません。そこで、階層ベイズを使うことでベータ分布のパラメータも推定するようにします。

階層ベイズによる比率の推定

事前分布をうまく設定することができれば、直感的によさそうな採用率を推定できることがわかりました。しかし、その事前分布のパラメータをどのように設定すればよいかわからないという問題があることもわかりました。ここではこの問題に対処するため、事前分布のパラメータにも事前分布を設定する、いわゆる階層ベイズを使うことで、採用率の事前分布のパラメータも推定します。

採用率の事前分布には今までと同じようにベータ分布を使いますが、ベータ分布のパラメータ\(\alpha\)、\(\beta\)にも事前分布を設定します。ここでは、後で扱いやすいように以下のような再パラメータ化を行います。\(\mu\)はベータ分布の期待値です。
\begin{eqnarray*}
\mu &=& \frac{\alpha}{\alpha + \beta} \\
\kappa &=& \alpha + \beta
\end{eqnarray*}

モデルとコードは以下のようになります。
f:id:livesense-analytics:20181004183644p:plain:w480

Rコード(クリックすると展開されます)

library(rstan)
library(tidyverse)

dat <- read_csv("sample.csv")

stan_dat <-
  list(M = nrow(dat),
       N = dat$N,
       R = dat$R)
 
fit <- stan(file = 'hier.stan',
            data = stan_dat,
            seed = 123)


Stanコード(hier.stan)(クリックすると展開されます)

data {
  int M;
  int<lower=0> R[M];  // 応募数(分母)
  int<lower=0> N[M];  // 採用数(分子)
}
 
parameters {
  vector<lower=0, upper=1>[M] theta;  // 採用率
  real<lower=0, upper=1> mu;          // 事前分布(ベータ分布)の期待値(= a / (a + b))
  real<lower=0> kappa;                // a + b
}
 
transformed parameters {
  real<lower=0> a;        // 事前分布(ベータ分布)のパラメータ
  real<lower=0> b;        // 事前分布(ベータ分布)のパラメータ
  a = mu * kappa;
  b = kappa * (1 - mu);
}
 
model {
  R ~ binomial(N, theta);
  theta ~ beta(a, b);     // 事前分布
}


このようにすることで、モデルとデータに合った、事前分布のパラメータ(の分布)が推定されます。ベータ分布のパラメータに固定値のような強い制約を課すのではなく、分布のような緩やかな制約を使っているところがポイントです。採用率の事前分布の条件が緩和されたことで、データに応じてパラメータが自動的に推定されるようになります。階層ベイズの場合も事前分布の事前分布を変更すると推定結果は影響を受けるものの、パラメータがモデルとデータから大きく逸脱した値をとる確率は低くなるため、階層化せず事前分布のパラメータ固定値を変更した時ほど大きな影響は受けることは少ないです。そのため、このようなケースで階層ベイズを使うと恣意的に推定結果を操作するのが難しくなるので、結果の客観性が高まるというメリットもあります。

推定結果は以下のようになります。

推定結果のサマリ(クリックすると展開されます)


前節と同様に推定結果を確認すると、階層ベイズでも直感的によさそうな採用率が推定されているように見えますね。前節で使ったベータ分布の固定パラメータは、この階層ベイズのパラメータを参考にして決めました。なお、ハイパーパラメータの有効サンプルサイズが小さいところは気になるので、必要に応じてサンプリング数を増やすなどの対応が必要です。

グループの違いを考慮した階層ベイズによる比率の推定

採用率は職種や地域によって大きく異なることが知られています。そこで最後に、職種や地域のようなグループごとの採用率の違いも考慮した推定をやってみましょう。

基本的には今までやってきたことと同じように、採用率の事前分布の事前分布のパラメータにさらに事前分布を設定します。グループごとに採用率を推定することもできますが、職種や地域によっては応募が少ないため全体傾向すら推定できないグループがあることが多いです。そこで、求人の採用率の事前分布にグループごとの採用率を使い、グループごとの採用率(の期待値)の事前分布にグループ別採用率の分布を設定します。このように階層を増やすことで、応募が少ないグループの採用率の推定にはグループ別採用率の分布が利用されるため、推定がしやすくなります。なお、ここではデータの分布の傾向については確認しませんが、グループごとの採用率分布も、応募の多いグループ内の求人の採用率分布も、前節までに見てきた全求人の採用率分布とほぼ同じような傾向になっているので、引き続き事前分布にはベータ分布を使います。

モデルとコードは以下のようになります。モデルが複雑になり計算が不安定になりやすいので、ハイパーパラメータ\(\kappa_g\)と\(\kappa_i\)に半正規分布を設定して安定化を図っています。また、Stan実行時にオプションadapt_deltaを0.99に設定しています。
f:id:livesense-analytics:20181004184637p:plain:w480

Rコード(クリックすると展開されます)

library(rstan)
library(tidyverse)

dat <- read_csv("sample.csv")

stan_dat <-
  list(M = nrow(dat),
       K = dat$G %>% unique %>% length,
       N = dat$N,
       R = dat$R,
       G = dat$G)

fit <- stan(file = 'hier2.stan', 
            data = stan_dat,
            seed = 123,
            control = list(adapt_delta = 0.99))


Stanコード(hier2.stan)(クリックすると展開されます)

data {
  int M;                            // データ数
  int K;                            // グループIDの最大値
  int<lower=0> N[M];                // 応募数(分母)
  int<lower=0> R[M];                // 採用数(分子)
  int<lower=1, upper=K> G[M];       // グループID
}

parameters {
  real<lower=0, upper=1> theta[M];  // 採用率
  
  vector<lower=0, upper=1>[K] mu;   // グループ別の採用率事前分布のパラメータ(期待値)
  vector<lower=0>[K] kappa;         // グループ別の採用率事前分布ののパラメータ(alpha + beta)
  
  real<lower=0, upper=1> mu_g;      // グループの事前分布のパラメータ(期待値)
  real<lower=0> kappa_g;            // グループの事前分布のパラメータ(alpha + beta)
}

transformed parameters {
  vector<lower=0>[K] a;             // グループ別の事前分布(ベータ分布)のパラメータ
  vector<lower=0>[K] b;             // グループ別の事前分布(ベータ分布)のパラメータ
  
  real<lower=0> a_g;                // グループの事前分布(ベータ分布)のパラメータ
  real<lower=0> b_g;                // グループの事前分布(ベータ分布)のパラメータ
  
  a = mu .* kappa;
  b = kappa .* (1 - mu);
  
  a_g = mu_g .* kappa_g;
  b_g = kappa_g .* (1 - mu_g);
}

model {
  R ~ binomial(N, theta);
  for (i in 1:M) {
    theta[i] ~ beta(a[G[i]], b[G[i]]);  // 採用率の事前分布(パラメータはグループ別)
  }
  mu ~ beta(a_g, b_g);                  // グループの事前分布(採用率の事前分布のパラメータの事前分布)
  kappa ~ normal(0, 100);
  kappa_g ~ normal(0, 100);
}


結果は以下のようになります。
推定結果のサマリ(クリックすると展開されます)


応募数と採用数が同じでも、グループによって推定採用率に違いがあることが確認できます。例として、推定採用率が比較的高いG=18に所属しているid=41の求人と、推定採用率が比較的低いG=24に所属しているid=35の求人について見てみましょう。どちらも応募1件採用0件ですが、推定採用率の期待値はそれぞれ0.215、0.1と異なっており、id=41のほうが高いことがわかります。期待通り、グループの違いが推定採用率に反映されていることがうかがえます。

一方で、データがほとんど変わらないのにパラメータが増えモデルが複雑になっているため、収束しにくくなっていることもわかると思います。ハイパーパラメータ\(\kappa_i\)の事前分布への依存性も高く、頑健性は高くありません。実務で利用する場合は、グループごとの違いが反映されるメリットと、計算が不安定化するデメリットを比較してどのモデルを使うかを決めることになると思います。また、より客観的にモデル選択をしたいのであれば、WBICを使う方法もあります。

最後に注意点を述べておこうと思います。今回のような推定結果を利用するときに気をつけなければいけいないことは、結果が事前分布や全体傾向に依存しているところです。事前分布を変更すると異なる推定値が得られることはよくありますし、データの抽出条件によって全体傾向が違うため抽出条件ごとに異なるモデルが必要とされることもあります。また、商品設計が変更されたり時間が経過したりしてデータの全体傾向が変わると、個別のデータの特徴には変化がなくとも事前分布が全体傾向に適合しなくなる可能性があるので注意が必要です。例えば、求人メディアの場合では、異なる採用率傾向をもつ求人が増えると全体傾向が変わる可能性があるので、事前分布の再検討が必要になるかもしれません。

まとめ

今回は階層ベイズを使って採用率などの分母が小さい比率を推定する方法を紹介しました。ポイントは以下の二点です

  1. 個別に見るとデータ量が少なくてまともな推定が難しそうでも、全体傾向などをうまく利用すると推定が可能になることがある
  2. 階層ベイズを利用することで事前分布のパラメータ(の分布)を自動的に推定できることがある

似たような方法は比率だけでなく別の小標本データにも応用できます。弊社にはベイズ統計に理解のあるディレクターもいるため、今回のようなベイズ推定結果を活用した取り組みも行われています。

なお、元ネタは統計モデリングの社内勉強会で階層ベイズを説明するときに使ったものです。この勉強会に参加したアナリストにもエンジニアにも受けがよかったのでブログ記事にしてみました。今後は小ネタでも役立ちそうなものがあれば紹介していこうと思います。