こんにちは、リブセンスで統計や機械学習関係の仕事をしている北原です。今回はレコメンデーションにも使える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}$（１サンプルのコンテキストは$\bf{x}$と表記）、モデルパラメータ$w_0 \in {\mathbb R}$、${\bf w} \in {\mathbb R}^n$、${\bf V} \in {\mathbb R}^{n \times k}$を使って、以下のようなモデルで表現されます。
モデル上のポイントは、交互作用項の係数$\hat{w}_{i, j}$がモデルパラメータ${\bf V}$の異なる行の内積になっているところにあります。これにより特徴量の組み合わせが学習データに存在しなくとも$\hat{w}_{i, j}$を計算できるため、レコメンデーションに利用しやすくなっているというのが１つ目の利点です。通常の回帰モデルでは、交互作用項の係数が特徴量の組み合わせの数$n(n+1)/2$だけ必要です。そのため、特徴量数$n$が多いとモデルパラメータ数が膨大になるという問題があります。しかし、FMでは${\bf V}$のサイズ分しか必要としません。通常は$k \ll n$となるよう因子数$k$を決めるので、同様の問題は起きにくくなるというのが２つ目の利点です。
一見すると交互作用項が複雑なのですが、各モデルパラメータについては線形になっているところが計算上のポイントです。二乗誤差を考えたとき、各モデルパラメータについては単純な二次関数になるため、簡単に最適値を計算することができます。
ここで後の説明をやりやすくするために、式($\ref{eq:fm}$)の別の表記を導入します。各モデルパラメータについて線形なので、着目するモデルパラメータを$\theta$としたとき式($\ref{eq:fm}$)は一般的に
と書くことができます。$g_{(\theta)}(\mathbf{x})$は$\theta$を含まない部分、$h_{(\theta)}(\mathbf{x})$は$\theta$にかかっている部分です。$\theta$以外のモデルパラメータはすでに与えられていて定数になっていると考えるとわかりやすいかもしれません。少し具体例を見てみましょう。例えば、特徴量数$n=3$、因子数$k=2$のとき、式($\ref{eq:fm}$)は
と書けます。$\theta=w_2$とすると
となります。同様に、$\theta=v_{3, 2}$とすると
となります。式($\ref{eq:lin}$)の表記を使うことで、$w$や$v$の違いにかかわらず統一的にモデルパラメータを扱えたり、微分を
のように簡潔に表記できます。しかし、これだけだと大変な計算が$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 descentやGradient descentに掲載されている図を見ると、直感的なイメージがつかめると思います。
SGDやGDのような単純な勾配法と比較したときのALSの利点は、学習率のような調整パラメータがないところです。近年では学習率を自動調整する方法がいろいろと提案されているので、あまり強いメリットを感じないかもしれません。それでも、最初から学習率を考えなくともよいのは利点の一つでしょう。
では、FMのモデルパラメータ更新式を導出しましょう。モデルパラメータ集合を$\Theta$、モデルパラメータ$\theta \in \Theta$の正則化係数を$\lambda_{(\theta)}$として正則化項を入れた二乗誤差
を最小化することでパラメータを推定します。モデルパラメータは複数ありますが、ALSでは１パラメータのみを変数として考え、残りを定数のように扱うところがポイントです。１つのモデルパラメータ$\theta$のみを変数としたとき、二乗誤差$L$を最小にする$\theta$の値は、$\theta$で微分したものを0とおいて整理すればよいですね。式($\ref{eq:lin}$)と式($\ref{eq:par_y}$)を使って
と変形し、$\theta$について整理すると
と、二乗誤差$L$を最小にする$\theta$を計算できます。式($\ref{eq:lin}$)を使っているのでシンプルな形にまとまっていますね。式($\ref{eq:fm}$)や式($\ref{eq:w}$)をそのまま愚直に微分してみると、違いが実感できると思います。
次に、１パラメータ１回の更新に必要な計算量について考えてみましょう。実際には全パラメータについて収束するまで値更新を繰り返すので、さらにパラメータ数と繰り返し数を乗じたものが全体の計算量になることに注意してください。
基本的には式($\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}$)の交互作用項は
と変形できるので、計算量は$O(s \, k \, n)$になります。さらに、データがスパースな場合、1サンプルあたりの非ゼロデータの個数の平均を$\overline{m}_R$とすると、$O(s \, k \, \overline{m}_R)$になります。レコメンデーションなどで使う場合は、$\overline{m}_R \ll n$となることがほとんどです。例えば、ユーザー数とアイテム数が数万のオーダーとして$n \sim O( 10^4 )$、１サンプルあたりの非ゼロデータ数が百程度として$\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$を使うと
と書けます。$e(\mathbf{x}_l, y_l | \Theta)$は$\theta$について線形なので、更新後のパラメータを$\theta^{\ast}$を含むパラメータ集合を$\Theta^{\ast}$とすると
と更新できます。更新前の$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}$)より
と計算できます。交互作用項の微分は式($\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$について線形なので
と更新でき、
と整理できます。つまり、$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$なので、計算の必要がありません。１パラメータに対応する非ゼロデータの個数の平均を$\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}$)を整理すると
となります。 これをコードに落とすと、モデルパラメータ推定のコードは以下のようになります。
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をレコメンデーションで使う場合は評価推定値計算の計算量も問題になりやすいです。なぜかというと、各ユーザーに対してレコメンド対象アイテムの評価推定値を計算するので、組み合わせが多くなってしまうためです。学習データで交差検証を行ったりコンペで使ったりしただけだと問題にならないので気づきにくいのですよね。この話題については、後日、別の記事で扱う予定です。
