Python: レコメンドの行列分解を確率的勾配降下法で実装する

過去3回の投稿で、行列分解(SVDとNMF)によるレコメンドを実装してきました。

ですが、いずれも欠測値やユーザごと・アイテムごとのバイアスに対応していませんでした。

今回の投稿では、欠測値やバイアスの概念を導入した行列分解について理解し、確率的勾配降下法で実装してみます。

欠測値とバイアスを考慮した行列分解(SVD)

欠測値とバイアスを考慮した行列分解について、Simon Funkのブログで投稿されています。元々はNetflix Prize(wikipedia)で考えられたアルゴリズムだそうです。

Netflix Update: Try This at Home

解説としてはこちらの論文が参考になります。

Y.Koren、R.Bell、C.Volinsky ： MATRIX FACTORIZATION TECHNIQUES FOR RECOMMENDER SYSTEM (2009) (PDF)
Y.Koren、R.Bell ： Advances in Collaborative Filtering (2011) (PDF)

Simonのアルゴリズムでは、ユーザuのアイテムiに対する評価値 $\hat{r_{ui}}$ は、以下式で予測します。

$\hat{r_{ui}}=\mu + b_i + b_u + q^T_i p_u$

μは全評価値の平均値
$b_i$ はアイテムiのバイアスを表しており、アイテムごとの高評価・低評価の差が反映されています
$b_u$ はユーザuのバイアスを表しており、ユーザごとの評価の甘さ・辛さの差が反映されています
qとuが分解された評価値行列となります
- qは、アイテムごとの特徴がk次元で表現された行列で、k行×M列(Mはアイテム数)
  - したがって、 $q_i$ は、アイテムiのk次元の特徴
- pは、ユーザごとの特徴がk次元で表現された行列で、k行×N列(Nはユーザ数)
  - したがって、 $p_u$ は、ユーザuのk次元の特徴

なお、推薦システムの界隈ではこの方法が一般的にSVDと呼ばれているようです。固有値から導出する数学的な特異値分解(SVD)と異なりますが、同じ"SVD"という用語が使われていますので、文脈で読み分ける必要があります。この投稿でも、SVDと呼称します。

確率的勾配降下法によるSVD

上の式では、バイアス( $b_i$ 、 $b_u$ )と分解後の行列(q, p)が未知のため、何らかの方法で導出する必要があります。

Simonの投稿では、確率的勾配降下法(SGD)が使われています。
SGDは、データ(ここでは評価値)をランダムに1つずつ取り出し、目的関数の値が小さくなるように変数を更新していく方法です。SGDについてはこちらの投稿が参考になります。

postd.cc

SGDで最小化を目指す目的関数は以下式となります。

$\min_{b, q, p} \sum_{(u, i) \in R} (r_{ui} - \mu - b_{i} - b_{u} - q^T_i p_u)^2 + \lambda (b_i^2 + b_u^2 + \mid q_i \mid ^2 + \mid p_u \mid ^2)$

第1項は、予測値と実際値の誤差( $r_{ui} - \hat{r_{ui}}$ )の2乗和です
第2項は、正則化項です
- ハイパーパラメータλは正則化制約の強さを表し、値が大きいほど強くなります
値がある評価値( $(u, i) \in R$ )のみで計算され、欠測値は除外されていることに注目してください

各値は以下の式で更新します。

$b_u \leftarrow b_u + \gamma (e_{ui} - \lambda b_u)$
$b_i \leftarrow b_i + \gamma (e_{ui} - \lambda b_i)$
$q_i \leftarrow q_i + \gamma (e_{ui} p_u - \lambda q_i)$
$p_u \leftarrow p_u + \gamma (e_{ui} q_i - \lambda p_u)$
$e_{ui}$ は誤差で、 $e_{ui} = r_{ui} - \mu - b_{i} - b_{u} - q^T_i p_u$
ハイパーパラメータγは学習率を表し、値が大きいほど1回の学習で値が大きく更新されます

実装

それではPythonで実装していきます。

例として、5本の映画に対して、4人のユーザが5段階(1〜5)で評価した、というケースを想定します(前回の投稿と同じです)。
未評価の値(欠測値)は0と置いておきます。この後の計算では除外されます。

f:id:ohke:20171223190013p:plain

import numpy as np
import pandas as pd
np.set_printoptions(precision=2, suppress=True)


# 評価値行列
x = np.matrix([[5, 5, 0, 3], [4, 4, 0, 0], [0, 3, 4, 0], [0, 1, 4, 5], [1, 0, 5, 4]])
# matrix([[5, 5, 0, 3],
#         [4, 4, 0, 0],
#         [0, 3, 4, 0],
#         [0, 1, 4, 5],
#         [1, 0, 5, 4]])

変数の初期化

q、p、 $b_u$ 、 $b_i$ を、-0.05〜+0.05の小さな乱数で初期化します(初期値は0でも良いかもしれません)。

def get_initial_values(x, k):
    # 各行列またはベクトルの値を、-0.05〜+0.05の乱数で初期化
    mu = np.sum(x) / np.count_nonzero(x)
    bu = (np.random.rand(x.shape[1]) - 0.5) * 0.1
    bi = (np.random.rand(x.shape[0]) - 0.5) * 0.1
    q = (np.random.rand(k, x.shape[0]) - 0.5) * 0.1
    p = (np.random.rand(k, x.shape[1]) - 0.5) * 0.1
    
    return (mu, bu, bi, q, p)


k = 2
get_initial_values(x, k)
# mu: 3.69
# bu:  array([-0.03,  0.  ,  0.  ,  0.04])
# bi:  array([-0.01, -0.02, -0.01, -0.05, -0.  ])
# q:
#  array([[-0.04, -0.01, -0.01,  0.02, -0.03],
#        [ 0.03, -0.04, -0.04, -0.02,  0.04]]),
# p:
#  array([[-0.04, -0.03,  0.  , -0.03],
#         [ 0.02,  0.02,  0.04,  0.02]]))

誤差の計算

誤差 $e_{ui}$ と目的関数の計算を実装します。

get_error_matrix関数は、1つの誤差の値が1要素となる行列にして返します
error_function関数は、目的関数の実装で、欠測値(0)はここで除外することで、学習へ影響を与えないようにします

def get_error_matrix(x, mu, q, p, bu, bi):
    # 誤差eの計算
    # 計算の簡略化のため、評価値0は除外していない
    return x - (mu + np.dot(q.T, p) + bu + np.matrix(bi).T) 


def error_function(x, mu, q, p, bu, bi, l):
    # 誤差に正則化制約を加えた目的関数を定義
    error_matrix = get_error_matrix(x, mu, q, p, bu, bi)
    error = np.sum(np.square(error_matrix[x > 0])) # 0を除外して誤差を計算
    
    regularization = l * (np.sum(np.square(bu)) + np.sum(np.square(bi)) + np.sum(np.square(q)) + np.sum(np.square(p)))
    
    return error + regularization


error_function(x, mu, q, p, bu, bi, 0.02)
# 23.03739901928402

SGDによる値の更新

最後にSGDで値を更新する実装です。

各エポックの最初に評価値をランダムに並び替えています(ここでも0は除外しています)。
そして、評価値を先頭から1つずつ取り出し、目的関数で値を求めて、各値を更新しています。

次元数kは2、エポック数は100の条件下で、目的関数の値は22.7から0.87まで下がります。

def sgd(x, k, epochs=100, l=0.02, gamma=0.02):
    # 確率的勾配法で分解した行列とバイアスを求める
    mu, bu, bi, q, p = get_initial_values(x, k)
    errors = [] # エポック毎の目的関数の値
    
    for epoch in range(epochs):
        error = error_function(x, mu, q, p, bu, bi, l)
        errors.append(error)
        
        # 欠測値を除いた評価値のインデックスを取得してランダムに並び替える
        x_i, x_u = np.where(x > 0)
        targets = np.arange(len(x_i))
        np.random.shuffle(targets)
        
        for target in targets:
            # 誤差行列を計算して、そのうち1つを選択
            error_matrix = get_error_matrix(x, mu, q, p, bu, bi)
            
            i = x_i[target]
            u = x_u[target]
            
            e_ui = error_matrix[i, u]
            
            # 分解した行列とバイアスを更新
            bu[u] = bu[u] + gamma * (e_ui - l * bu[u])
            bi[i] = bi[i] + gamma * (e_ui - l * bi[i])
            q[:, i] = q[:, i] + gamma * (e_ui * p[:, u] - l * q[:, i])
            p[:, u] = p[:, u] + gamma * (e_ui * q[:, i] - l * p[:, u])

    error = error_function(x, mu, q, p, bu, bi, l)
    errors.append(error)
    print('Error: {}'.format(error))
    # Error: 0.8715459314064704

    expected = mu + bu + np.matrix(bi).T + np.dot(q.T, p)
    
    return expected, errors
    

expected, errors = sgd(x, 2)

得られた値は以下になります。アイテムやユーザの特徴・嗜好というよりは、アイテムやユーザごとのバイアスからの影響が強そうです。

$b_u$ では、ユーザA(1列目)のユーザが辛め、ユーザC(3列目)のユーザが甘めの傾向が反映されています
$b_i$ では、タイタニック(2行目)のアイテムが高評価、サマーウォーズ(4行目)のアイテムが低評価となっています
qは、1列目を見ると、ローマの休日(1行目)とサマーウォーズ(4行目)、AKIRA(5行目)が近くなっているようです
pは、1行目を見るとユーザA(1列目)とB(2列目)、また、2行目ではBとC(3列目)が近いことを示しているようです

0で置き換えていた欠測値にも、予測された評価値が入っており、例えばユーザAとBのアニメーション映画に対する評価は低い、といった嗜好が反映されています。

$$ \hat{R}=\mu + b_u + b_i + q^Tp = 3.7+\left( \begin{matrix} -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \end{matrix} \right)+\left( \begin{matrix} 0.7 & 0.7 & 0.7 & 0.7 \\ 0.8 & 0.8 & 0.8 & 0.8 \\ -0.4 & -0.4 & -0.4 & -0.4 \\ -0.7 & -0.7 & -0.7 & -0.7 \\ -0.6 & -0.6 & -0.6 & -0.6 \end{matrix} \right)+\left( \begin{matrix} 1.2 & 0.4 \\ -0.3 & 0.0 \\ -0.2 & 0.0 \\ 1.1 & 0.7 \\ 0.9 & 0.2 \end{matrix} \right) \left( \begin{matrix} -1.2 & -0.9 & 0.4 & 1.1 \\ 0.0 & 0.5 & 0.2 & -0.6 \end{matrix} \right)= \left( \begin{matrix} 4.9 & 5.1 & 5.4 & 3.1 \\ 3.9 & 4.0 & 5.6 & 4.6 \\ 2.7 & 2.8 & 4.3 & 3.3 \\ 1.2 & 1.1 & 4.1 & 4.9 \\ 1.4 & 2.0 & 4.5 & 4.1 \end{matrix} \right) $$

エポック毎に目的関数の値をプロットしました。エポック毎に0に近づくことがわかります。

f:id:ohke:20171229215509p:plain