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

前々回の投稿( Python: レコメンドの行列分解を確率的勾配降下法で実装する - け日記 )では、欠測値やバイアスを考慮した行列分解を、確率的勾配降下法(SGD)で求める実装を行いました。

今回はおなじみMovieLens100K Datasetへ、このアルゴリズムを適用します。

性能の問題

前々回の投稿のおさらいです。

ohke.hateblo.jp

SGDでは以下の目的関数を最小化する4つのパラメータ( $b_u$ , $b_i$ , q, p)を求めていました。

$\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つずつに対して、誤差の計算とパラメータの更新(下式)を行っていました。

$b_i \leftarrow b_i + \gamma (e_{ui} - \lambda b_i)$
$b_u \leftarrow b_u + \gamma (e_{ui} - \lambda b_u)$
$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}$ も再計算しているためです。評価値が100,000個あれば誤差も100,000個あるため、1回の再計算のコストも大変高くなります。

実際、前々回の実装をそのままMovieLens100K Datasetに適用したところ、学習に時間がかかりすぎてしまいました。

ミニバッチの導入

そこで、ミニバッチを導入します。

ミニバッチ勾配降下法(MSGD、単にSGDとも)では、評価値1つずつに対して誤差の計算とパラメータの更新を行うのではなく、例えば100個ごとで誤差の計算とパラメータの更新をまとめて行います。これにより、誤差 $e_{ui}$ の計算回数が1エポックあたり100,000回から1,000回まで削減できます。
各種勾配降下法の比較については、こちらの投稿が参考になります。

勾配降下法の最適化アルゴリズムを概観する | コンピュータサイエンス | POSTD

確率的勾配降下法でパラメータを学習する関数sgd()を、ミニバッチを導入して書き直したmsgd()を作成します。それ以外の関数の実装は、前々回と同じです。

ミニバッチで一度に扱う評価値の個数を表すbatch_sizeを引数に追加
delta_bu, delta_bi, delta_q, delta_pにミニバッチ内での更新値が累積される

import numpy as np
import pandas as pd
import math
import time
import codecs
from scipy import sparse
from matplotlib import pyplot


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)


def get_error_matrix(x, mu, q, p, bu, bi):
    # 誤差eの計算
    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]))
    
    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


def msgd(x, k, epochs=100, l=0.02, gamma=0.02, batch_size=1):
    # 確率的勾配法で分解した行列とバイアスを求める
    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)
        rating_count = len(x_i)
        batch_count = math.ceil(rating_count / batch_size)
        targets = np.arange(rating_count)
        np.random.shuffle(targets)
        
        # ミニバッチ単位で誤差の計算とパラメータの更新を行う
        for batch_base in range(batch_count):
            error_matrix = get_error_matrix(x, mu, q, p, bu, bi)
            
            delta_bu = np.zeros_like(bu)
            delta_bi = np.zeros_like(bi)
            delta_q = np.zeros_like(q)
            delta_p = np.zeros_like(p)
            
            for batch_offset in range(batch_size):
                # 対象の評価値のインデックスを計算
                target_index = batch_size * batch_base + batch_offset
                i = x_i[target_index]
                u = x_u[target_index]
                e_ui = error_matrix[i, u]
                
                # パラメータの更新値を累積
                delta_bu[u] += gamma * (e_ui - l * bu[u])
                delta_bi[i] += gamma * (e_ui - l * bi[i])
                delta_q[:, i] += gamma * (e_ui * p[:, u] - l * q[:, i])
                delta_p[:, u] += gamma * (e_ui * q[:, i] - l * p[:, u])
                
            # パラメータの更新
            bu += delta_bu
            bi += delta_bi
            q += delta_q
            p += delta_p

    error = error_function(x, mu, q, p, bu, bi, l)
    print('Error: {}'.format(error))
    
    # 予測した評価値を生成
    expected = mu + bu + np.matrix(bi).T + np.dot(q.T, p)
    
    return expected, errors

MovieLens100K Datasetでの評価

今回もMovieLens 100K Datasetを使います。

以下からダウンロードして、今回作成するPythonファイルと同じディレクトリに解凍しておきます。以前の投稿( word2vecでitem2vecを実装して映画を推薦する - け日記 )も参考にしてみてください。

MovieLens | GroupLens

評価値はnumpy行列にして持ちます。欠測値は0で補完しています。アイテムIDとタイトルの一覧もDataFrameにロードしておきます。
また、指定したユーザに対して未評価のアイテムの評価値を予測してランキングを返すpredict_ranking()と、それを表示するprint_ranking()を実装しています。

# 評価値を持つファイルはu.data。楽に行列化するために、まずはDataFrameでロードする
ratings = pd.read_csv('ml-100k/u.data', delimiter='\t', header=None).iloc[:, :3]
ratings.rename(columns={0: 'user_id', 1: 'item_id', 2: 'rating'}, inplace=True)
print('ratings.shape: {}'.format(ratings.shape)) # ratings.shape: (100000, 3)

# item_idを行、user_idを列とする1682×943の行列(欠測値は0埋め)
rating_matrix = np.asmatrix(ratings.pivot(index='item_id', columns='user_id', values='rating').fillna(0), dtype=np.float64)
print('rating_matrix.shape: {}'.format(rating_matrix.shape)) # rating_matrix.shape: (1682, 943)

# pd.read_csvで読み込もうとすると文字コードエラーになるため、
# codecs.openでエラーを無視してファイルを読み込んだ後にDataFrameにする
with codecs.open('ml-100k/u.item', 'r', 'utf-8', errors='ignore') as f:
    item_df = pd.read_table(f, delimiter='|', header=None).ix[:, 0:1]
    item_df.rename(columns={0: 'item_id', 1: 'item_title'}, inplace=True)
    item_df.set_index('item_id', inplace=True)

items = pd.Series(item_df['item_title'])
print('items.shape: {}'.format(items.shape)) # items.shape: (1682, 1)


def predict_ranking(user_index, reducted_matrix, original_matrix, n):
    # 対象ユーザの予測評価値
    reducted_vector = np.array(reducted_matrix[:, user_index].T)[0]
    
    # 評価済みのアイテムの値は0にする
    filter_vector = np.array(original_matrix[:, user_index].T)[0] == 0
    predicted_vector = np.multiply(reducted_vector, filter_vector)

    # 上位n個のアイテムのインデックスを返す
    return [(i, predicted_vector[i]) for i in np.argsort(predicted_vector)[::-1][:n]]


def print_ranking(user_ids, items, reducted_matrix):
    for user_id in user_ids:
        print('User ID: {}'.format(user_id))
        predicted_ranking = predict_ranking(user_id - 1, reducted_matrix, rating_matrix, 10)
        for item_id, rating in predicted_ranking:
            # アイテムID, 映画タイトル, 予測した評価値を表示
            print('{}: {} [{}]'.format(item_id + 1, items[item_id + 1], rating))

それでは予測評価値の行列を作ります。

k(qとpの行サイズ)は30とし、これまでのSVDやNMFと同じ値に設定
エポック数は1000として、エポック毎に目的関数の値(誤差)をプロットすることで、十分低減できているかを確認
- 学習に約90分を要しました

# 1000エポックにして誤差の傾向を確認
%time reducted_matrix, errors = msgd(rating_matrix, k=30, epochs=1000, batch_size=1000, gamma=0.001)
# Error: 22387.638436133657
# CPU times: user 1h 22min 20s, sys: 19min 38s, total: 1h 41min 59s
# Wall time: 1h 29min 37s

エポック毎の目的関数の値をプロットしたグラフが以下です。
エポック数600を超えたあたりから大きな変化はないため、エポック数1000で十分と判断できます。最終的な誤差は22387となり、k=30ではこのくらいが限界と予想されます。

f:id:ohke:20180108132848p:plain

最後に、2人のユーザへ推薦させてみます。

SVDやNMFとの結果と比較するために、前回と同じ2人のユーザ(ユーザID 267と868)を対象にランキングを生成しました。この2人は"Akira(1988)"と"Ghost in the Shell(Kokaku kidotai)(1995)"を5で評価しています(つまりSF+アニメーション好きです)。

akira_item_id = 206
ghost_in_the_shell_item_id = 1240

# Akira (1988)とGhost in the Shell (Kokaku kidotai) (1995)を5で評価したユーザIDを取得
user_ids = np.where(np.multiply(rating_matrix[akira_item_id - 1] == 5, rating_matrix[ghost_in_the_shell_item_id - 1] == 5))[1] + 1

print_ranking(user_ids, items, reducted_matrix)

結果を見ると、以下のように考察できるかと思います。

SVDとNMFと同じタイトルがいくつか見られます("Clerks (1994)"や"One Flew Over the Cuckoo's Nest (1975)")
- ですが、SVDやNMFと比較すると、SF要素が弱く、納得感が低いように見えます(バイアスの影響が大きかったのかもしれません)
ユーザ同士で同じタイトルが推薦されています("Paths of Glory (1957)", "Clerks (1994)", Third Man, The (1949)など)
- 類似したユーザと捉えることはできている
267の方が全体的に予測した評価値が高くなっています
- ユーザのバイアスの影響を受けています(267の平均評価値が3.9に対して、868は3.0)

User ID: 267
42: Clerks (1994) [5.343020693286839]
516: Local Hero (1983) [5.306247713944661]
1449: Pather Panchali (1955) [5.267915793013931]
511: Lawrence of Arabia (1962) [5.217895974837784]
357: One Flew Over the Cuckoo's Nest (1975) [5.210585560183143]
960: Naked (1993) [5.165327044088528]
641: Paths of Glory (1957) [5.1401291180844675]
1512: World of Apu, The (Apur Sansar) (1959) [5.123277190111157]
513: Third Man, The (1949) [5.109739675120826]
652: Rosencrantz and Guildenstern Are Dead (1990) [5.06174853656031]

User ID: 868
641: Paths of Glory (1957) [6.020237492858888]
522: Down by Law (1986) [5.515328289793725]
171: Delicatessen (1991) [5.411841155332155]
513: Third Man, The (1949) [5.308531517941752]
42: Clerks (1994) [5.307793583401524]
262: In the Company of Men (1997) [5.251078640118084]
430: Duck Soup (1933) [5.223552176891947]
836: Ninotchka (1939) [5.205202414859825]
443: Birds, The (1963) [5.204741054839308]
519: Treasure of the Sierra Madre, The (1948) [5.189719602700117]

まとめ

欠測値やバイアスを考慮した行列分解(SVD)によるレコメンドを実装し、MovieLens100Kに適用しました。
また、確率的勾配降下法にミニバッチを導入することで、巨大な評価値行列を扱う実用ケースにおいても、なんとか耐えうる実装にしました。