ここでは，行列分解 (matrix factorization) による協調フィルタリングを実装してみます．まず，Pandas，NumPy に加えて scipy.optimize ライブラリをインポートします．さらに，Numpyでの小数点以下の表示桁数を調整するとともに，指数表示にならないように出力形式を設定しておきます．
import pandas as pd
import numpy as np
import scipy.optimize as optimize
np.set_printoptions(precision=3) # 小数点以下の表示桁数を設定
np.set_printoptions(suppress=True) # 指数表示を行わないように
次に，GitHub のリポジトリから行列形式のレイティングデータをデータフレームに読み込みます．
url = "https://github.com/rinsaka/sample-data-sets/blob/master/collaborative_filtering_matrix.csv?raw=true";
df = pd.read_csv(url)
print(df)
id name A B C D E 0 0 Eto 5.0 4.0 3.0 NaN NaN 1 1 Sato NaN 2.0 4.0 2.0 3.0 2 2 Kato 2.0 3.0 4.0 NaN NaN 3 3 Muto 2.0 3.0 2.0 NaN 2.0 4 4 Kito 4.0 3.0 3.0 4.0 NaN 5 5 Goto 2.0 5.0 NaN 2.0 3.0 6 6 Bito 3.0 3.0 2.0 NaN 1.0 7 7 Saito 4.0 4.0 5.0 NaN 5.0 8 8 Naito 5.0 5.0 4.0 2.0 1.0 9 9 Koto 3.0 NaN 2.0 2.0 2.0
目的は顧客 Kato（ と Eto ）に対して商品「D」と商品「E」のどちらを優先的におすすめすべきかを決めることです．顧客の好みが異なれば自ずとおすすめする商品も異なるはずです．
ここでは，行列分解と呼ばれる協調フィルタリングを考えます．この手法ではおおむね次のようなステップで顧客に対するおすすめ商品を決定します．
まず顧客によるレイティング行列 \(\bf{Y}\) が与えられ，\(i\) 行 \(j\) 列要素である \(Y_{ij}\) が顧客 \(i\) が商品 \(j\) のレイティングであるとします．ただし，レイティング行列の要素のすべてが観測されている必要はありません．
顧客 \(i\) が商品 \(j\) に対して与えるスコアを \(s_{ij}\) として，次のように仮定します． \begin{eqnarray} s_{ij} &=& {\boldsymbol{u}}_i^{\prime} {\boldsymbol{v}}_j \end{eqnarray} ここで，\({\boldsymbol{u}}_i\) は顧客 \(i\) の潜在因子 (latent factor) で，\({\boldsymbol{v}}_j\) は商品 \(j\) の潜在因子です．また，\({\boldsymbol{u}}_i\)，\({\boldsymbol{v}}_j\) はいずれも \(L\) 次元のベクトルです．
行列因子分解では次の最適化問題を解くことで \({\boldsymbol{u}}_i\) と \({\boldsymbol{v}}_j\) の最尤推定ができるようになります． \begin{eqnarray} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \argmin_{{\boldsymbol{u}}_i, {\boldsymbol{v}}_j, \forall_i \forall_j} \sum_{(i, j) \in \Omega} \left( y_{ij} - {\boldsymbol{u}}_i^{\prime} {\boldsymbol{v}}_j \right)^2 \end{eqnarray}
なお，上の方法では過学習してしまう傾向があります．したがって，多くの場合 \(L_2\) ノルムによる正則化が用いられます． \begin{eqnarray} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \argmin_{{\boldsymbol{u}}_i, {\boldsymbol{v}}_j, \forall_i \forall_j} \sum_{(i, j) \in \Omega} \left( y_{ij} - {\boldsymbol{u}}_i^{\prime} {\boldsymbol{v}}_j \right)^2 + \lambda_1 \sum_{i} || {\boldsymbol{u}}_i ||^2 + \lambda_2 \sum_{j} || {\boldsymbol{v}}_j ||^2 \end{eqnarray}
これから行列因子分解の手順を確認しながらコードを記述していきます．
すでに読み込んでいる Pandas のデータフレームからレイティングに関する列だけを取得して NumPy 配列に格納し，レイティング行列 y を生成します．
y = df.iloc[:, 2:].values
print(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.]]
顧客数 \(M\) は行列 \(\bf Y\) の行数，アイテム数 \(N\) は列数に対応するので，これらの情報を取得しておきます．
M, N = y.shape
print("顧客数 M =", M)
print("アイテム数 N =",N)
顧客数 M = 10 アイテム数 N = 5
次に潜在因子数 \(L\) をとりあえず 3 に設定します．この値は最終的には十分に議論する必要があるでしょう．
L = 3
print("潜在因子数 L =", L)
潜在因子数 L = 3
まず，顧客に関する仮の潜在因子行列を1で初期化し，顧客 2 の潜在因子ベクトル \({\boldsymbol{u}}_2\) だけ別の値を設定します．
U = np.ones((M, L))
U[2, 0] = 1
U[2, 1] = 2
U[2, 2] = 3
print(U)
[[1. 1. 1.]
[1. 1. 1.]
[1. 2. 3.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]
同様に商品に関する仮の潜在因子行列を 1 で初期化し，商品ID=0 （商品A）と商品ID=1 （商品B）の潜在因子ベクトル \({\boldsymbol{v}}_0\), \({\boldsymbol{v}}_1\) についてだけ別の値を設定します．
V = np.ones((N, L))
V[0, 0] = 3; V[0, 1] = 1; V[0, 2] = -1
V[1, 0] = -2; V[1, 1] = 1; V[1, 2] = 3
print(V)
[[ 3. 1. -1.] [-2. 1. 3.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.]]
2つの潜在因子行列からスコアは \begin{eqnarray} s_{ij} &=& {\boldsymbol{u}}_i^{\prime} {\boldsymbol{v}}_j \end{eqnarray} を使って計算できます．ここで，\({\boldsymbol{u}}_2 = \left(\begin{array}{c}1&2&3\end{array} \right)\) は行列 \(\bf U\) の \(2\) 行目の要素である（ただし先頭が0行目としています）ので，\({\boldsymbol{u}}_2^{\prime}\) はこれを転置した \begin{eqnarray} {\boldsymbol{u}}_2^{\prime} = \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right) \end{eqnarray} となります．また，\({\boldsymbol{v}}_0 = \left(\begin{array}{c}3&1&-1\end{array} \right)\) は行列 \(\bf V\) の \(0\) 行目の要素です．したがって，\(s_{2,0}\) は \begin{eqnarray} s_{2,0} &=& {\boldsymbol{u}}_2^{\prime} {\boldsymbol{v}}_0 \\ &=& \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right) \left(\begin{array}{c}3&1&-1\end{array} \right) \\ &=& 1 \times 3 + 2 \times 1 + 3 \times (-1) \\ &=& 2 \end{eqnarray} となります．つまり，顧客2の商品「A」のスコアは 2 となりました．さらに，\(s_{2,1}\) は \begin{eqnarray} s_{2,1} &=& {\boldsymbol{u}}_2^{\prime} {\boldsymbol{v}}_1 \\ &=& \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right) \left(\begin{array}{c}-2&1&3\end{array} \right) \\ &=& 1 \times (-2) + 2 \times 1 + 3 \times 3 \\ &=& 9 \end{eqnarray} となります．
これは，商品ごとに実際には観測できない3種類の属性（特性・特徴など）があり，商品「A」の属性が 3, 1, -1 で商品「B」の属性が -2, 1, 3 であることを意味しています．また，顧客にも実際には観測できない3種類の属性に対する反応の重みがあり，顧客2の Kato さんは 1, 2, 3 の重みで商品「A」や「B」を評価することを意味しています．この結果，Kato さんの商品「A」に対するスコアは 2, 商品「B」に対するスコアは 9 となったことを意味しています．
Python では次のように導出することができます．
s = np.dot(U, V.T)
print(s)
[[3. 2. 3. 3. 3.]
[3. 2. 3. 3. 3.]
[2. 9. 6. 6. 6.]
[3. 2. 3. 3. 3.]
[3. 2. 3. 3. 3.]
[3. 2. 3. 3. 3.]
[3. 2. 3. 3. 3.]
[3. 2. 3. 3. 3.]
[3. 2. 3. 3. 3.]
[3. 2. 3. 3. 3.]]
レイティングとその推定値の差は次にように計算できます．
print(y-s)
[[ 2. 2. 0. nan nan]
[nan 0. 1. -1. 0.]
[ 0. -6. -2. nan nan]
[-1. 1. -1. nan -1.]
[ 1. 1. 0. 1. nan]
[-1. 3. nan -1. 0.]
[ 0. 1. -1. nan -2.]
[ 1. 2. 2. nan 2.]
[ 2. 3. 1. -1. -2.]
[ 0. nan -1. -1. -1.]]
2乗誤差を求めてみます．
print((y-s)**2)
[[ 4. 4. 0. nan nan] [nan 0. 1. 1. 0.] [ 0. 36. 4. nan nan] [ 1. 1. 1. nan 1.] [ 1. 1. 0. 1. nan] [ 1. 9. nan 1. 0.] [ 0. 1. 1. nan 4.] [ 1. 4. 4. nan 4.] [ 4. 9. 1. 1. 4.] [ 0. nan 1. 1. 1.]]
2乗誤差について行ごとの，つまり，顧客ごとの和を計算します．
print(np.nansum((y-s)**2, axis=1))
[ 8. 2. 40. 4. 3. 11. 6. 13. 19. 3.]
2乗誤差について列ごとの，つまり，アイテムごとの和を計算します．
print(np.nansum((y-s)**2, axis=0))
[12. 65. 13. 5. 14.]
2乗誤差について総和，つまり残差平方和を計算します．これは 109.0 となりました．
print(np.nansum((y-s)**2))
109.0
これは目的関数である \begin{eqnarray} \sum_{(i, j) \in \Omega} \left( y_{ij} - {\boldsymbol{u}}_i^{\prime} {\boldsymbol{v}}_j \right)^2 \end{eqnarray} を求めたことになります．つまり，この2乗誤差の総和が最小化されるように \({\boldsymbol{u}}_i\) と \({\boldsymbol{v}}_j\) を決めれば良いということになります．
上で求めた2乗誤差の総和が最小化されるように最適な \({\boldsymbol{u}}_i\) と \({\boldsymbol{v}}_j\) を求めるために，scipy.optimize.minimize() を利用します．このとき，scipy.optimize.minimize() に与える独立変数は1次元配列の形式でなければなりません．今回は2つの行列を与えたいので，この2つの行列を連結した上で1次元配列に変形します．
まず行列 \({\bf{U}}\) の形状を確認すると，10行3列であることがわかります．
print(U)
[[1. 1. 1.] [1. 1. 1.] [1. 2. 3.] [1. 1. 1.] [1. 1. 1.] [1. 1. 1.] [1. 1. 1.] [1. 1. 1.] [1. 1. 1.] [1. 1. 1.]]
次に，行列 \({\bf{V}}\) の形状は5行3列です．
print(V)
[[ 3. 1. -1.] [-2. 1. 3.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.]]
行列 \({\bf{U}}\) と 行列 \({\bf{V}}\) の列数が等しいのでこれらを連結します．
print(np.concatenate([U,V]))
[[ 1. 1. 1.] [ 1. 1. 1.] [ 1. 2. 3.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.] [ 3. 1. -1.] [-2. 1. 3.] [ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.]]
さらに，これを1次元配列に変形して UV に代入します．これで45個の要素からなる1次元配列となりました．
UV = np.concatenate([U,V]).ravel()
print(UV)
[ 1. 1. 1. 1. 1. 1. 1. 2. 3. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3. 1. -1. -2. 1. 3. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
すでに目的関数を求める手順は示しましたが，これを関数化します．このとき，独立変数は1次元配列として受け取るようにして，関数内部で2つの行列に戻すようにしておきます．
def get_rss_UV(UV, Y, M, N, L):
"""
目的関数
行列Uと行列Vを一つの一次元配列にした UV を受け取り，行列 U と V に分割・変換して演算する
"""
u, v = np.split(UV, [M*L])
u = u.reshape(M, L)
v = v.reshape(N, L)
S = np.dot(u, v.T)
return np.nansum((Y-S)**2)
いま定義した関数を呼び出して残差平方和 (residual sum of squares) が正しく計算できることを確認しておきます．この結果，すでに求めた 109.0 と同じ結果が得られました．
get_rss_UV(UV, y, M, N, L)
109.0
では，残差平方和が最小となるような最適化を行います．まず，行列 \({\bf{U}}\), \({\bf{V}}\) の初期値は \([-1, 1)\) の一様乱数で与えることにします．
# 乱数列の初期化
# rng = np.random.default_rng()
rng = np.random.default_rng(seed=1)
# -1, 1 の一様乱数
U = 2.0 * rng.random((M, L)) - 1.0
V = 2.0 * rng.random((N, L)) - 1.0
UV = np.concatenate((U.ravel(), V.ravel()))
print(UV)
[ 0.024 0.901 -0.712 0.897 -0.376 -0.153 0.655 -0.182 0.099 -0.945 0.507 0.076 -0.341 0.577 -0.394 -0.093 -0.732 -0.194 -0.593 -0.475 0.501 -0.439 -0.03 0.961 0.923 0.45 0.082 -0.446 -0.679 0.94 0.032 -0.768 0.247 0.553 0.226 0.835 -0.921 0.057 -0.081 -0.875 0.283 0.705 0.186 -0.48 0.68 ]
なお，初期値での目的関数値は 431.76 です．
get_rss_UV(UV, y, M, N, L)
431.76386648930963
scipy.optimize.minimize() を使って最適化します．このとき，最適化アルゴリズムにはネルダーミード法を使い，許容誤差 (tol) は 0.0001，最大繰り返し回数 (maxiter) は500万回に設定しています．これらは適宜調整が必要であると思われます．この結果，最適化が成功し，目的関数の最小値が 12.2078 となりました．
results = optimize.minimize(
get_rss_UV, UV, args=(y, M, N, L),
method='Nelder-Mead',
tol=0.0001,
options={'maxiter': 5000000}
)
# print(results) # すべての結果を表示する場合はここを有効に
print(' fun:', results.fun)
print('message:', results.message)
print(' status:', results.status)
print('success:', results.success)
fun: 12.207832863842437 message: Optimization terminated successfully. status: 0 success: True
得られた解とスコア，並びに RSS（残差平方和），MSE（Mean Squared Error：平均二乗誤差） を表示します．この結果，スコア \(s\) はレイティング行列 \(y\) に概ね近い値になっていることがわかります．しかしながら，レイティング行列 \(y\) が空値であった箇所については絶対値が極端に大きな値になっています．これは学習データに過剰に適合してしまっている過学習が起こっていることを意味します．この原因は，行列 \({\bf{U}}\), \({\bf{V}}\) の最適解の一部に絶対値の大きな値があることです．
u, v = np.split(results.x, [M*L])
u = u.reshape(M, L)
v = v.reshape(N, L)
print('----- u -------')
print(u)
print('----- v -------')
print(v)
s = np.dot(u, v.T)
print('----- s -------')
print(s)
print('----- y -------')
print(y)
print('----- RSS -------')
print(np.nansum((y-s)**2))
print('----- RSS -------')
print(np.nansum((y-s)**2))
----- u ------- [[ -1.845 15.255 2.03 ] [ 10.708 7.24 7.158] [ 69.442 -15.984 6.157] [ 6.418 6.806 -3.556] [ 16.571 5.537 4.637] [ -0.631 10.982 11.798] [ 1.544 10.036 -10.915] [ 17.608 1.385 74.815] [ 40.186 1.65 -0.337] [ -9.368 11.162 12.387]] ----- v ------- [[ 0.101 0.303 0.022] [ 0.106 0.295 0.024] [ 0.103 0.221 0.04 ] [-0.237 5.899 -5.334] [ 0.028 0.199 0.057]] ----- s ------- [[ 4.474 4.349 3.259 79.598 3.102] [ 3.432 3.442 2.981 1.991 2.154] [ 2.325 2.823 3.841 -143.598 -0.858] [ 2.63 2.606 2.019 57.596 1.336] [ 3.454 3.505 3.107 3.998 1.838] [ 3.52 3.448 2.829 2.003 2.84 ] [ 2.952 2.866 1.939 117.055 1.423] [ 3.852 4.041 5.09 -395.038 5.029] [ 4.558 4.755 4.473 2.003 1.453] [ 2.703 2.585 1.996 2.001 2.661]] ----- 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 ------- 12.207832863842437 ----- MSE ------- 0.3130213554831394
ここでは，正則化と呼ばれる手法によって，過学習を抑える方法について説明します．ここでは \(L_2\) ノルムによる正則化項を加えた目的関数を考えます． \begin{eqnarray} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \argmin_{{\boldsymbol{u}}_i, {\boldsymbol{v}}_j, \forall_i \forall_j} \sum_{(i, j) \in \Omega} \left( y_{ij} - {\boldsymbol{u}}_i^{\prime} {\boldsymbol{v}}_j \right)^2 + \lambda_1 \sum_{i} || {\boldsymbol{u}}_i ||^2 + \lambda_2 \sum_{j} || {\boldsymbol{v}}_j ||^2 \end{eqnarray} これにより \({\bf{U}}\), \({\bf{V}}\) の推定値の絶対値も小さくするように働きます．
目的関数を次のように定義します．
def get_L2rss_UV(UV, Y, M, N, L, lam1, lam2):
"""
L2正則化項を加えた目的関数
UとVを一つの一次元配列にした UV を受け取り，行列 U と V に変換して演算する
"""
u, v = np.split(UV, [M*L])
u = u.reshape(M, L)
v = v.reshape(N, L)
S = np.dot(u, v.T)
# L2 ノルム
ul2 = np.linalg.norm(u, ord=2)
vl2 = np.linalg.norm(v, ord=2)
return np.nansum((Y-S)**2 + lam1 * ul2 + lam2 * vl2)
当初設定した \({\bf{U}}\), \({\bf{V}}\) の値を使って新たな目的関数の値を計算してみます．なお，\(\lambda_1\) と \(\lambda_2\) には 0.01 (=1.0E-02)という値を設定しています．これも今後調整が必要と考えられるパラメータです．具体的にはクロスバリデーションなどの方法が利用して，\(\lambda_1\) と \(\lambda_2\) に加えて \(L\) についても調整すると良いでしょう．
U = np.ones((M, L))
U[2, 0] = 1; U[2, 1] = 2; U[2, 2] = 3
V = np.ones((N, L))
V[0, 0] = 3; V[0, 1] = 1; V[0, 2] = -1
UV = np.concatenate((U.ravel(), V.ravel()))
lam1 = 1.0E-02
lam2 = 1.0E-02
get_L2rss_UV(UV, y, M, N, L, lam1, lam2)
68.08783560041476
実際に最適化を行ってみます．この結果，\({\bf{U}}\), \({\bf{V}}\) の推定値の絶対値に極端に大きな値がなくなり，スコアにも極端な値がなくなったように感じられます．さらに，RSS や MSE の値も改善していることがわかります．なお，次のページでは \(\lambda_1\) と \(\lambda_2\) に 0.001 (=1.0E-03) という値を使って更によい解が得られたことを示しています．
# 乱数列の初期化
# rng = np.random.default_rng()
rng = np.random.default_rng(seed=1)
# -1, 1 の一様乱数
U = 2.0 * rng.random((M, L)) - 1.0
V = 2.0 * rng.random((N, L)) - 1.0
UV = np.concatenate((U.ravel(), V.ravel()))
lam1 = 1.0E-02
lam2 = 1.0E-02
results = optimize.minimize(
get_L2rss_UV, UV, args=(y, M, N, L, lam1, lam2),
method='Nelder-Mead',
tol=0.001,
options={'maxiter': 5000000}
)
print(' fun:', results.fun)
print('message:', results.message)
print(' status:', results.status)
print('success:', results.success)
u, v = np.split(results.x, [M*L])
u = u.reshape(M, L)
v = v.reshape(N, L)
print('----- u -------')
print(u)
print('----- v -------')
print(v)
s = np.dot(u, v.T)
print('----- s -------')
print(s)
print('----- y -------')
print(y)
print('----- RSS -------')
print(np.nansum((y-s)**2))
print('----- MSE -------')
print(np.nansum((y-s)**2)/np.count_nonzero(~np.isnan(y)))
fun: 13.406992209328342 message: Optimization terminated successfully. status: 0 success: True ----- u ------- [[ 0.041 -4.811 0.681] [-2.231 -1.086 0.576] [-1.819 -1.249 0.546] [-1.451 -0.13 0.589] [-2.699 -0.926 1.056] [-0.597 5.405 1.925] [-1.507 -0.999 0.826] [-4.347 -0.042 1.384] [-2.137 -2.009 1.48 ] [-1.655 -0.278 1.012]] ----- v ------- [[ 0.044 -0.639 2.873] [ 0.14 -0.307 3.463] [-1.27 -0.702 -0.251] [-1.626 0.436 -0.554] [-1.142 0.413 0.081]] ----- s ------- [[ 5.034 3.842 3.154 -2.543 -1.98 ] [ 2.25 2.016 3.45 2.833 2.144] [ 2.287 2.02 3.049 2.11 1.605] [ 1.712 1.877 1.786 1.976 1.651] [ 3.508 3.565 3.811 3.397 2.784] [ 2.049 4.921 -3.519 2.262 3.071] [ 2.946 2.957 2.407 1.556 1.374] [ 3.812 4.199 5.201 6.281 5.057] [ 5.444 5.446 3.752 1.777 1.73 ] [ 3.014 3.36 2.042 2.008 1.856]] ----- 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 ------- 8.391952104006844 ----- MSE ------- 0.21517825907709856
なお，このページでは小さなデータセットについて，\({\bf{U}}\) と \({\bf{V}}\) を一気に最適化する方法を採用しました．しかしこれは大きなデータセットではうまく動作しない可能性があります．実際上は交互最小二乗法や確率的勾配降下法などの最適化手法を利用する必要があるかもしれません．また，このままのコードでは計算時のメモリ消費量が非常に大きく，計算自体に失敗することになるでしょう．これらの詳細は参考図書を参照してください．一方で，Surprise ライブラリを利用することで，大きなデータであってもより簡単に顧客間類似度のモデルや行列分解のモデルを実行することができるようになります．