以前のページでは，行列分解の仕組みを理解しながら，そのモデル化と最適化の Python プログラムを実装しました．しかしながら，実装したアルゴリズムは効率的な書き方になっていない部分もあることから，大きなデータで実行すると膨大なメモリを消費するとともに，長時間の計算が必要になることでしょう．ここでは，Surprise という推薦システムのライブラリの使用例を示します．
まず，pip list コマンドなどで surprise がインストールされていることを確認してください．インストールがまだの場合は，前のページの手順でインストールしてください．
Surprise のインストールができれば，インポートします．前のページとの違いは KNNWithMeans ではなく SVD をインポートしていることだけです．
import pandas as pd
import numpy as np
from surprise import Dataset, Reader, SVD
np.set_printoptions(precision=3) # 小数点以下の表示桁数を設定
np.set_printoptions(suppress=True) # 指数表示を行わないように
この後の作業はしばらく前のページと同様です．サンプルデータを GitHub から読み込み，データフレームに格納します．このとき，行列形式のデータではなく，次のような形式のデータを読み込んでいることに注意指定ください．
url_ratings = "https://github.com/rinsaka/sample-data-sets/blob/master/collaborative_filtering_ratings.csv?raw=true"
df_ratings = pd.read_csv(url_ratings)
print(df_ratings)
customer item rating 0 8 3 2 1 1 1 2 2 6 1 3 3 9 2 2 ...（中略）... 35 7 0 4 36 9 0 3 37 4 3 4 38 4 1 3
surprise ではレイティングのスケールをあらかじめ設定する必要があります．Reader を使って rating の最小値 1 と最大値 5を次のように指定します．ここで設定したスケールは後のステップにおいてレイティングの予測値を求めたときの下限・上限としても利用されることに注意してください．
reader = Reader(rating_scale=(1, 5))
列名を指定してデータフレームからデータセットを読み込みます．
data = Dataset.load_from_df(df_ratings[["customer", "item", "rating"]], reader)
読み込んだ data にどのような形式で値が格納されているか確認しておきます．
vars(data)
{'reader': <surprise.reader.Reader at 0x11785b210>,
'has_been_split': False,
'df': customer item rating
0 8 3 2
1 1 1 2
2 6 1 3
3 9 2 2
...（中略）...
35 7 0 4
36 9 0 3
37 4 3 4
38 4 1 3,
'raw_ratings': [(8, 3, 2.0, None),
(1, 1, 2.0, None),
(6, 1, 3.0, None),
(9, 2, 2.0, None),
...（中略）...
(7, 0, 4.0, None),
(9, 0, 3.0, None),
(4, 3, 4.0, None),
(4, 1, 3.0, None)]}
今回はすべてのレイティングデータを学習データとして設定します．
trainset = data.build_full_trainset()
設定した学習データがどのような形式になっているかも確認しておきます．
vars(trainset)
{'ur': defaultdict(list,
{0: [(0, 2.0), (3, 5.0), (4, 1.0), (1, 5.0), (2, 4.0)],
1: [(1, 2.0), (2, 4.0), (0, 2.0), (4, 3.0)],
2: [(1, 3.0), (3, 3.0), (4, 1.0), (2, 2.0)],
3: [(2, 2.0), (0, 2.0), (4, 2.0), (3, 3.0)],
4: [(2, 4.0), (1, 3.0), (3, 2.0)],
5: [(2, 3.0), (3, 4.0), (0, 4.0), (1, 3.0)],
6: [(2, 2.0), (1, 3.0), (4, 2.0), (3, 2.0)],
7: [(3, 5.0), (1, 4.0), (2, 3.0)],
8: [(4, 3.0), (0, 2.0), (3, 2.0), (1, 5.0)],
9: [(2, 5.0), (1, 4.0), (4, 5.0), (3, 4.0)]}),
'ir': defaultdict(list,
{0: [(0, 2.0), (8, 2.0), (1, 2.0), (3, 2.0), (5, 4.0)],
1: [(1, 2.0),
(2, 3.0),
(4, 3.0),
(7, 4.0),
(9, 4.0),
(8, 5.0),
(6, 3.0),
(0, 5.0),
(5, 3.0)],
2: [(3, 2.0),
(4, 4.0),
(1, 4.0),
(5, 3.0),
(6, 2.0),
(9, 5.0),
(7, 3.0),
(0, 4.0),
(2, 2.0)],
3: [(2, 3.0),
(0, 5.0),
(7, 5.0),
(5, 4.0),
(8, 2.0),
(4, 2.0),
(6, 2.0),
(9, 4.0),
(3, 3.0)],
4: [(8, 3.0),
(0, 1.0),
(2, 1.0),
(6, 2.0),
(1, 3.0),
(3, 2.0),
(9, 5.0)]}),
'n_users': 10,
'n_items': 5,
'n_ratings': 39,
'rating_scale': (1, 5),
'_raw2inner_id_users': {8: 0,
1: 1,
6: 2,
9: 3,
2: 4,
4: 5,
3: 6,
0: 7,
5: 8,
7: 9},
'_raw2inner_id_items': {3: 0, 1: 1, 2: 2, 0: 3, 4: 4},
'_global_mean': None,
'_inner2raw_id_users': None,
'_inner2raw_id_items': None}
上の結果の中に，学習データ内のユーザ数やアイテム数が格納されていることがわかったので，その値を表示してみます．
print(f"number of trainset customers = {trainset.n_users}")
print(f"number of trainset items = {trainset.n_items}")
number of trainset customers = 10 number of trainset items = 5
また，学習データは元のデータと必ずしも同じ並び順ではありません．次のコードは必ずしも必要ありませんが，学習データを行列形式に変換したものと，元の順番に並べたものを表示してみます．学習データと元のデータの並び順の関係は，上の結果にある _raw2inner_id_users や _raw2inner_id_items を手がかりに得られます．
trainset_matrix = np.zeros((trainset.n_users, trainset.n_items))
raw_matrix = np.zeros((trainset.n_users, trainset.n_items)) # 元のデータ
trainset_matrix[:,:] = np.nan
raw_matrix[:,:] = np.nan
for u, vals in trainset.ur.items():
for i, rating in vals:
# print(u, i, rating)
trainset_matrix[u, i] = rating
raw_matrix[trainset.to_raw_uid(u), trainset.to_raw_iid(i)] = rating
print(trainset_matrix)
print("-"*25)
print(raw_matrix)
[[ 2. 5. 4. 5. 1.] [ 2. 2. 4. nan 3.] [nan 3. 2. 3. 1.] [ 2. nan 2. 3. 2.] [nan 3. 4. 2. nan] [ 4. 3. 3. 4. nan] [nan 3. 2. 2. 2.] [nan 4. 3. 5. nan] [ 2. 5. nan 2. 3.] [nan 4. 5. 4. 5.]] ------------------------- [[ 5. 4. 3. nan nan] [nan 2. 4. 2. 3.] [ 2. 3. 4. nan nan] [ 2. 3. 2. nan 2.] [ 4. 3. 3. 4. nan] [ 2. 5. nan 2. 3.] [ 3. 3. 2. nan 1.] [ 4. 4. 5. nan 5.] [ 5. 5. 4. 2. 1.] [ 3. nan 2. 2. 2.]]
ここまでの手順は前のページとほぼ同じでした．ここから，行列分解を行います．具体的には，（潜在）因子数，学習のエポック数，学習率などを設定し，SVD （特異値分解 :
Singular Value Decomposition の意味）でモデルを設定し，fit で当てはめを行います．なお，use_biase = False とすることで，バイアスを利用しないモデルを利用します．これを use_biase = True にするとバイアスを利用するモデルが利用されるので，さらに良い当てはめ結果になることが予想されます（これも是非試してください）．また，正則化パラメータの意味はこちらで確認してください．
# 因子数
n_factors = 3
# エポック数
n_epochs = 100000
# バイアスを利用するか
biased = False
# 学習率
lr_all = 0.005
# 正則化パラメータ
reg_all = 0.02
# モデルの設定
mf = SVD(n_factors=n_factors,
n_epochs=n_epochs,
biased=biased,
lr_all=lr_all,
reg_all=reg_all
)
# モデルの当てはめ（学習）を行う
mf.fit(trainset)
<surprise.prediction_algorithms.matrix_factorization.SVD at 0x10641e2d0>
学習ができれば，顧客（ユーザ）IDと商品（アイテム）IDを指定することで predict モジュールによってレイティング予測値を求めることができます．例えば Kato さん (uid=2) の商品 D (iid=3) についてのレイティング予測値を求めます．この結果から予測値は est = 2.38 であることが読み取れます．
uid = 2
iid = 3
pred = mf.predict(uid, iid, verbose=False)
print(pred)
user: 2 item: 3 r_ui = None est = 2.38 {'was_impossible': False}
特定の顧客と商品の組み合わせについては求められたので，すべての顧客と商品の組み合わせについて，レイティングの予測値を求めて行列形式のデータ形式にします．さらに，予測値 (s)，学習データ (y) と残差平方和 (RSS)，平均2乗誤差 (MSE) を表示します．この結果，精度良く予測ができていることがわかり，RSS や MSE も以前のページの結果よりもさらに小さな値になっていることがわかります．
predict_matrix = np.zeros((trainset.n_users, trainset.n_items))
for uid in range(trainset.n_users):
for iid in range(trainset.n_items):
pred = mf.predict(uid, iid, verbose=False)
predict_matrix[uid, iid] = pred.est
print("----- s -------")
print(predict_matrix)
print("----- y -------")
print(raw_matrix)
print("----- RSS -------")
print(np.nansum((predict_matrix-raw_matrix)**2))
print("----- MSE -------")
print(np.nansum((predict_matrix-raw_matrix)**2)/np.count_nonzero(~np.isnan(raw_matrix)))
----- s ------- [[4.936 4.016 2.981 3.1 2.521] [1. 2.111 3.852 1.919 3.067] [1.972 3.048 3.927 2.38 2.998] [2.26 2.452 2.425 1.754 1.784] [3.992 2.991 3.01 3.901 3.978] [2.056 4.873 5. 2.059 2.933] [3.017 2.873 2.074 1.43 1. ] [3.922 4.097 4.899 4.376 4.97 ] [4.852 5. 3.801 1.938 1.176] [2.923 2.419 2.009 2.114 1.902]] ----- y ------- [[ 5. 4. 3. nan nan] [nan 2. 4. 2. 3.] [ 2. 3. 4. nan nan] [ 2. 3. 2. nan 2.] [ 4. 3. 3. 4. nan] [ 2. 5. nan 2. 3.] [ 3. 3. 2. nan 1.] [ 4. 4. 5. nan 5.] [ 5. 5. 4. 2. 1.] [ 3. nan 2. 2. 2.]] ----- RSS ------- 0.8645121743948669 ----- MSE ------- 0.02216697883063766
行列分解による予測ができるようになりましたが，ここで予測値がどの様に計算されているかを理解します．具体的な数式はこちらで確認できますが，モデルの当てはめを mf.fit で行なった後に mf クラスにどのような変数が格納されているかを確認します．この中にある pu と qi が行列分解によって分解された行列です．
vars(mf)
{'n_factors': 3,
'n_epochs': 100000,
'biased': False,
'init_mean': 0,
'init_std_dev': 0.1,
'lr_bu': 0.005,
'lr_bi': 0.005,
'lr_pu': 0.005,
'lr_qi': 0.005,
'reg_bu': 0.02,
'reg_bi': 0.02,
'reg_pu': 0.02,
'reg_qi': 0.02,
'random_state': None,
'verbose': False,
'bsl_options': {},
'sim_options': {'user_based': True},
'trainset': <surprise.trainset.Trainset at 0x1178670d0>,
'bu': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
'bi': array([0., 0., 0., 0., 0.]),
'pu': array([[-0.488, 2.028, -1.502],
[-1.759, -0.398, -0.662],
[-0.378, 1.266, -0.627],
[-0.789, 1.149, -0.09 ],
[-1.582, 0.421, -0.69 ],
[-1.717, 1.435, 0.463],
[-0.848, 0.79 , -0.427],
[-0.98 , 2.03 , -0.276],
[-1.793, 0.285, -2.044],
[-2.377, 1.161, -0.147]]),
'qi': array([[-1.431, 0.887, 0.372],
[-0.999, 1.316, -1.324],
[-1.767, 0.428, -1.379],
[-0.606, 2.114, -0.179],
[-1.941, 0.348, 0.316]])}
まず，顧客（ユーザ）に関する潜在変数を意味する行列を表示します．
print(mf.pu)
[[-0.488 2.028 -1.502] [-1.759 -0.398 -0.662] [-0.378 1.266 -0.627] [-0.789 1.149 -0.09 ] [-1.582 0.421 -0.69 ] [-1.717 1.435 0.463] [-0.848 0.79 -0.427] [-0.98 2.03 -0.276] [-1.793 0.285 -2.044] [-2.377 1.161 -0.147]]
次に商品（アイテム）に関する行列を表示します．
print(mf.qi)
[[-1.431 0.887 0.372] [-0.999 1.316 -1.324] [-1.767 0.428 -1.379] [-0.606 2.114 -0.179] [-1.941 0.348 0.316]]
この2つの行列について積を計算します．しかしながら，これは surprise 内部の順序で並んでいるため，直感的に理解することは困難です．なお，@ は行列の積を求めるための演算子で，.T は行列の転置を行うメソッドです．
s_pred = mf.pu @ mf.qi.T
print(s_pred)
[[1.938 5.145 3.801 4.852 1.176] [1.919 2.111 3.852 0.344 3.067] [1.43 2.873 2.074 3.017 0.975] [2.114 2.419 2.009 2.923 1.902] [2.38 3.048 3.927 1.972 2.998] [3.901 2.991 3.01 3.992 3.978] [1.754 2.452 2.425 2.26 1.784] [3.1 4.016 2.981 4.936 2.521] [2.059 4.873 6.109 2.056 2.933] [4.376 4.097 4.899 3.922 4.97 ]]
より理解がしやすいように元のデータと同じ順番に並べ替えた行列 s_raw_pred を作成します．この上で，行列の積で求めた s_raw_pred，predict で得られた結果，学習データ，行列の積と predict の差を表示します．
s_raw_pred = np.zeros((s_pred.shape))
for u, row in enumerate(s_pred):
for i, rating in enumerate(row):
# print(u, i, rating)
s_raw_pred[trainset.to_raw_uid(u), trainset.to_raw_iid(i)] = rating
print("----- matrix product -------")
print(s_raw_pred)
print("----- predict -------")
print(predict_matrix)
print("----- y -------")
print(raw_matrix)
print("----- difference -------")
print(predict_matrix - s_raw_pred)
----- matrix product ------- [[4.936 4.016 2.981 3.1 2.521] [0.344 2.111 3.852 1.919 3.067] [1.972 3.048 3.927 2.38 2.998] [2.26 2.452 2.425 1.754 1.784] [3.992 2.991 3.01 3.901 3.978] [2.056 4.873 6.109 2.059 2.933] [3.017 2.873 2.074 1.43 0.975] [3.922 4.097 4.899 4.376 4.97 ] [4.852 5.145 3.801 1.938 1.176] [2.923 2.419 2.009 2.114 1.902]] ----- predict ------- [[4.936 4.016 2.981 3.1 2.521] [1. 2.111 3.852 1.919 3.067] [1.972 3.048 3.927 2.38 2.998] [2.26 2.452 2.425 1.754 1.784] [3.992 2.991 3.01 3.901 3.978] [2.056 4.873 5. 2.059 2.933] [3.017 2.873 2.074 1.43 1. ] [3.922 4.097 4.899 4.376 4.97 ] [4.852 5. 3.801 1.938 1.176] [2.923 2.419 2.009 2.114 1.902]] ----- y ------- [[ 5. 4. 3. nan nan] [nan 2. 4. 2. 3.] [ 2. 3. 4. nan nan] [ 2. 3. 2. nan 2.] [ 4. 3. 3. 4. nan] [ 2. 5. nan 2. 3.] [ 3. 3. 2. nan 1.] [ 4. 4. 5. nan 5.] [ 5. 5. 4. 2. 1.] [ 3. nan 2. 2. 2.]] ----- difference ------- [[ 0. 0. 0. 0. 0. ] [ 0.656 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. ] [ 0. 0. -1.109 0. 0. ] [ 0. 0. 0. 0. 0.025] [ 0. 0. 0. 0. 0. ] [ 0. -0.145 0. 0. 0. ] [ 0. 0. 0. 0. 0. ]]
上の結果を確認すると，最後の「difference」から「matrix product」と「predict」において，4つの要素で値が異なっていることが読み取れます．これらの値の違いを注意深く観察すると，「matrix product」の結果のうち，上で設定した1から5というレイティングの範囲に収まっていない値が，下限の1や上限の5に変更されていることがわかります．
したがって，予測値を範囲内に収めるための関数 check_rating_scale を定義し，s_raw_pred 行列を作成するときにこの関数を呼び出すように変更します．
def check_rating_scale(rating, min_scale=1.0, max_scale=5.0):
"""
評価予測値の範囲を制限する
"""
if rating < min_scale:
return min_scale
if rating > max_scale:
return max_scale
return rating
s_raw_pred = np.zeros((s_pred.shape))
for u, row in enumerate(s_pred):
for i, rating in enumerate(row):
# print(u, i, rating)
s_raw_pred[trainset.to_raw_uid(u), trainset.to_raw_iid(i)] = check_rating_scale(rating)
print("----- matrix product -------")
print(s_raw_pred)
print("----- predict -------")
print(predict_matrix)
print("----- y -------")
print(raw_matrix)
print("----- difference -------")
print(predict_matrix - s_raw_pred)
----- matrix product ------- [[4.936 4.016 2.981 3.1 2.521] [1. 2.111 3.852 1.919 3.067] [1.972 3.048 3.927 2.38 2.998] [2.26 2.452 2.425 1.754 1.784] [3.992 2.991 3.01 3.901 3.978] [2.056 4.873 5. 2.059 2.933] [3.017 2.873 2.074 1.43 1. ] [3.922 4.097 4.899 4.376 4.97 ] [4.852 5. 3.801 1.938 1.176] [2.923 2.419 2.009 2.114 1.902]] ----- predict ------- [[4.936 4.016 2.981 3.1 2.521] [1. 2.111 3.852 1.919 3.067] [1.972 3.048 3.927 2.38 2.998] [2.26 2.452 2.425 1.754 1.784] [3.992 2.991 3.01 3.901 3.978] [2.056 4.873 5. 2.059 2.933] [3.017 2.873 2.074 1.43 1. ] [3.922 4.097 4.899 4.376 4.97 ] [4.852 5. 3.801 1.938 1.176] [2.923 2.419 2.009 2.114 1.902]] ----- y ------- [[ 5. 4. 3. nan nan] [nan 2. 4. 2. 3.] [ 2. 3. 4. nan nan] [ 2. 3. 2. nan 2.] [ 4. 3. 3. 4. nan] [ 2. 5. nan 2. 3.] [ 3. 3. 2. nan 1.] [ 4. 4. 5. nan 5.] [ 5. 5. 4. 2. 1.] [ 3. nan 2. 2. 2.]] ----- difference ------- [[0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.]]
これにより分解された行列の積によって求めたレイティングの予測値と surprise の predict メソッドを利用して求めたレイティングの予測値が一致していることが確認できたので，安心して predict メソッドを利用できそうです．なお，ここでは，use_biase = False として当てはめを実行したので，mf.qi.T と mf.qi.T からレイティングの予測値を計算できました．一方で，use_biase = True として当てはめた場合は，次の通り，mf.bu， mf.bi 及び，trainset._global_mean も利用しなければならないことに注意してください（それでも範囲内に収まらない場合は右辺全体に check_rating_scale を併用するとよいでしょう） ．具体的な数式はやはりこのページを参照してください．
s_raw_pred = np.zeros((s_pred.shape))
for u, row in enumerate(s_pred):
for i, rating in enumerate(row):
# s_raw_pred[trainset.to_raw_uid(u), trainset.to_raw_iid(i)] = check_rating_scale(rating)
s_raw_pred[trainset.to_raw_uid(u), trainset.to_raw_iid(i)] = \
trainset._global_mean + mf.bu[u] + mf.bi[i] + rating