行列分解によるレコメンドを整理する #2~Matorix Factorizationによるレコメンド~

0. はじめに

　この記事は、機械学習によるレコメンド手法の基礎について、私が理解しやすいと考えた順序・考え方で整理し直すものになります。対象としては「データ分析や機械学習の基本手法を理解しているものの、レコメンドに関しては未経験である」という人を想定しています。

　今回は、機械学習の手法を用いた行列分解であるMatrix Factorizationについて整理します。

0-0. 特異値分解による行列分解レコメンドの簡単な説明

　本記事では、行列分解によるレコメンドの基本的な概念や機械学習の基礎（勾配降下法など）については既知のものとして扱います。行列分解によるレコメンドについては前回の記事でまとめたので、よろしければご覧ください。

　
　簡単に説明すると、「評価値行列」を「ユーザ因子行列」と「アイテム因子行列」に分解し、行列の積や行列自体を利用する手法を「行列分解」と呼びます。

その1つの手法として「特異値分解」によるレコメンドがあるのですが、この手法には、

• 未観測の値（NaN）を0などで補完する必要があり、未観測の評価値を予測するのには使えない
• アイテム間の序列付けには使えるが、未観測を0などで補完するため、バイアスがある

という欠点があります。これらの欠点を解消しているのが、今回紹介するMatrix Factorizationという手法です。

1. 行列分解による評価値予測を一般化

　前回同様、以下のような評価値行列を考えます。

user_id/item_id	1	2	3	4	5
1	0.0	4.0	4.5	NaN	NaN
2	3.5	2.0	5.0	2.5	5.0
3	2.5	4.0	0.5	2.5	4.0
4	1.0	3.5	2.5	3.5	2.5
5	4.0	4.5	0.5	1.5	3.0
6	3.0	4.5	3.5	2.0	1.0

　特異値分解による行列分解は、例えば以下のようになります。

特異値分解では、「評価値行列を分解してユーザ・アイテム因子行列を求める」という発想で、評価値行列を解析的に（数学的な手法によって）分解して上記のような計算結果を得ました。

　この結果は、ユーザ因子$p$、アイテム因子$q$、評価値$r$を用いて、以下のように一般化できます。

$r$においては、1つ目のインデックスがユーザID、2つ目のインデックスがアイテムIDを表しています。$p$、$q$においては1つ目のインデックスがユーザ・アイテムのID、2つ目のインデックスは因子IDを表しています。例えば、$r_{11}$は「ユーザ1のアイテム1に対する評価値」です。この評価値は、ユーザ1の因子ベクトル$p_{1} = (p_{11} \ p_{12} \ p_{13} \ p_{14} \ p_{15})$（ユーザ因子行列の1行目）とアイテム1の因子ベクトル$q_{1} = (q_{11} \ q_{12} \ q_{13} \ q_{14} \ q_{15})^T$（アイテム因子行列の1列目、$T$は転置）の積によって、以下のように表すことができます。

r_{11} = p_{11} \times q_{11} + p_{12} \times q_{12} + p_{13} \times q_{13} + p_{14} \times q_{14} + p_{15} \times q_{15}

2. Matorix Factorization

　Matrix Factorizationの理論を整理します。この手法は、ユーザ・アイテム因子行列を、教師あり学習によって導出するという手法です。手法の元論文としてはよくYehuda Koren, Robert Bell and Chris Volinsky(2009)が参照されています。

2-1. 評価値を正解データにした学習

　まずは、線形回帰における教師あり学習について認識を合わせます。最もシンプルな定式化を考えると、正解データ$t$、予測値$y$、入力$x$、パラメータ$w$、バイアス$b$を用いて、以下のように表せると思います。

\begin{align}
&回帰式：y = b + w_1x_1 + w_2x_2 + … \\
&誤差関数：E = \frac{1}{2}(t - y)^2 \\
\end{align}

この式におけるパラメータ$w_j$（j = 1, 2, ...）をSGD（確率的勾配降下法）で最適化する場合、更新式は学習率$η$を用いて以下になります。

\begin{align}
&w_2の更新式：w_1 ← w_1 - η\frac{∂E}{∂w_1} \\
&w_2の更新式：w_2 ← w_2 - η\frac{∂E}{∂w_2} \\
&...
\end{align}

　この定式化を、前セクションで示したユーザ・アイテム因子行列の積に適用します。
　あるユーザuのアイテムiに対する評価値を$r_{ui}$とし、ユーザ・アイテム因子行列（ベクトル）の積による予測値を便宜上$\hat{r_{uj}}$と置きます。上記の線形回帰の例になぞらえると、以下のように定式化できます。

\begin{align}
&ユーザ1, アイテム1の回帰式：\hat{r_{11}} = p_{11} \times q_{11} + p_{12} \times q_{12} + p_{13} \times q_{13} + p_{14} \times q_{14} + p_{15} \times q_{15}\\
&ユーザ1, アイテム1の誤差関数：E_{11} = \frac{1}{2}(r_{11} - \hat{r_{11}})^2 \\
\\
&ユーザ1, アイテム2の回帰式：\hat{r_{12}} = p_{11} \times q_{21} + p_{12} \times q_{22} + p_{13} \times q_{23} + p_{14} \times q_{24} + p_{15} \times q_{25}\\
&ユーザ1, アイテム2の誤差関数：E_{12} = \frac{1}{2}(r_{12} - \hat{r_{12}})^2 \\
&...\\
\\
&ユーザu, アイテムiの回帰式：\hat{r_{ui}} = p_{u1} \times q_{i1} + p_{u2} \times q_{i2} + p_{u3} \times q_{i3} + p_{u4} \times q_{i4} + p_{u5} \times q_{i5}\\
&ユーザu, アイテムiの誤差関数：E_{ui} = \frac{1}{2}(r_{ui} - \hat{r_{ui}})^2 \\
&...\\
\\
\end{align}

$E_{ui}$は、ユーザuのアイテムiにおける誤差関数です。このように、ユーザ数×アイテム数分の回帰式が定義されます。

　次に、更新式について考えます。式に$x$（入力）がないため混乱するかもしれませんが、一度どちらかの因子（例えばアイテム因子$q$）を$x$と見なして式を眺めると分かりやすいかもしれません（SGDによる更新を行う場合は、$p$、$q$どちらにもランダムな初期値が格納されています）。
　例えばユーザ1のユーザ因子$p_{1k}$（k=1,2,3,4,5）の更新式は以下のようになります。

\begin{align}
&p_{11}の更新式：\\
&p_{11} ← p_{11} - η\frac{∂E_{11}}{∂p_{11}} \\
&p_{11} ← p_{11} - η\frac{∂E_{12}}{∂p_{11}} \\
&...\\
&p_{11} ← p_{11} - η\frac{∂E_{15}}{∂p_{11}} \\
\\
&p_{12}の更新式：\\
&p_{12} ← p_{12} - η\frac{∂E_{11}}{∂p_{12}} \\
&p_{12} ← p_{12} - η\frac{∂E_{12}}{∂p_{12}} \\
&...\\
&p_{12} ← p_{12} - η\frac{∂E_{15}}{∂p_{12}} \\
\\
&...\\
\\
&p_{15}の更新式：\\
&p_{15} ← p_{15} - η\frac{∂E_{11}}{∂p_{15}} \\
&p_{15} ← p_{15} - η\frac{∂E_{12}}{∂p_{15}} \\
&...\\
&p_{15} ← p_{15} - η\frac{∂E_{15}}{∂p_{15}} \\
\end{align}

各ユーザ因子$p_{11}, p_{12}, ...$の更新式がアイテムの数だけ存在します。すべての回帰式を書き出してみればわかると思いますが、ユーザ1のユーザ因子は「ユーザ1のアイテム1に関する回帰式」から「ユーザ1のアイテム5に関する回帰式」まですべてに登場します。つまり、各ユーザ因子$p$は、1回のエポックでアイテム数分更新されるということになります。
　同様に、アイテム1のアイテム因子$q_{ik}$（k=1,2,3,4,5）の更新式は以下のようになります。

\begin{align}
&q_{11}の更新式：\\
&q_{11} ← q_{11} - η\frac{∂E_{11}}{∂q_{11}} \\
&q_{11} ← q_{11} - η\frac{∂E_{21}}{∂q_{11}} \\
&...\\
&q_{11} ← q_{11} - η\frac{∂E_{61}}{∂q_{11}} \\
\\
&q_{12}の更新式：\\
&q_{12} ← q_{12} - η\frac{∂E_{11}}{∂q_{12}} \\
&q_{12} ← q_{12} - η\frac{∂E_{21}}{∂q_{12}} \\
&...\\
&q_{12} ← q_{12} - η\frac{∂E_{61}}{∂q_{12}} \\
\\
&...\\
\\
&q_{15}の更新式：\\
&q_{15} ← q_{15} - η\frac{∂E_{11}}{∂q_{15}} \\
&q_{15} ← q_{15} - η\frac{∂E_{21}}{∂q_{15}} \\
&...\\
&q_{15} ← q_{15} - η\frac{∂E_{61}}{∂q_{15}} \\
\end{align}

各アイテム因子$q$は、1回のエポックでユーザ数分更新されていることが分かります。

　このように、各因子（パラメータ）は、1回のエポックでも複数回更新されていますが、SGDの場合、当然1回目のエポックの更新値は初期値に依存しています。そのため、数エポック更新を繰り返すことで、$p$、$q$を徐々に最適化していきます。

2-2. 未観測な評価値に対する予測

　前セクションのような計算をしていくと、「ユーザ1のアイテム4（$r_{14}=NaN$）」や「ユーザ1のアイテム5（$r_{11}=NaN$）」に対する更新式はどのように定義すればよいのか？という疑問にあたります。当然、正解データが欠損していては勾配降下法は使えません。だからと言って、欠損値を補完してしまっては、特異値分解の時と同じ課題が生じてしまいます。

　結論としては、欠損値の補完は不要であり、またこれら未観測の評価値はパラメータの更新に利用する必要はありません。なぜなら、「ユーザ1のアイテム4」や「ユーザ1のアイテム5」における更新をせずとも、目的である「ユーザ・アイテム因子行列の導出」はすでに完了しているからです。

　例えば、$r_{14}$を予測するためにはユーザ1の因子ベクトル$p_{1} = (p_{11} \ p_{12} \ p_{13} \ p_{14} \ p_{15})$とアイテム4の因子ベクトル$q_{1} = (q_{11} \ q_{12} \ q_{13} \ q_{14} \ q_{15})^T$が必要になります。今回の評価値行列では、ユーザ1はアイテム1、2、3に対して評価値があるため、$r_{11}$、 $r_{12}$、$r_{13}$を正解データとして用いることで$p_{1}$を最適化できます。同様に、アイテム4はユーザ2、3、4、5からの評価値があるため、$r_{42}$、 $r_{43}$、$r_{44}$、$r_{45}$を正解データとして用いることで$q_{4}$を最適化できます。

　以下は、実際にMatrix Factorizationを実施した例です。

試しにユーザ1のアイテム2の予測評価値を計算してみると

\begin{align}
& \hat{r_{12}} = -1.687 \times -2.186 + -0.304 \times -0.034 + 0.222 \times 1.375 \\
&　　\,\,\, + -0.129 \times -0.312 + -0.129 \times -0.476 \\
&　\,\,\, = 4.1050 \\
&　\,\,\, ≒ 4.0
\end{align}

$r_{12}=4.0$なので、かなり良く予測できているように思われます。
　
　また、すでにお気づきかとは思いますが、ユーザ・アイテム因子行列が導出できているということは、評価値が欠損しているユーザ1のアイテム4・5の評価値を予測することもできます。

\begin{align}
& \hat{r_{14}} = -1.687 \times -1.048 + -0.304 \times -0.182 + 0.222 \times 1.062 \\
&　　\,\,\, + -0.129 \times 0.617 + -0.129 \times 0.452 \\
&　\,\,\, = 1.9212 \\
&　\,\,\, ≒ 2.0 \\
\\
& \hat{r_{15}} = -1.687 \times -0.352 + -0.304 \times -1.449 + 0.222 \times 1.610 \\
&　　\,\,\, + 0.467 \times 0.617 + -0.129 \times 0.026 \\
&　\,\,\, = 1.3281 \\
&　\,\,\, ≒ 1.5 \\
\end{align}

今回のMatrix Factorizationに基づくと「ユーザ1は、アイテム4に2.0の、アイテム5に1.5の評価値を付けると予測される」という結果になりました。

　以上より、Matrix Factorizationによる行列分解を行うことで、欠損値を補完することなくユーザ・アイテム行列を求めることができ、さらにその欠損値（未観測の評価値）を予測するタスクにも対応できるということが理解できたと思います。

2-3. 因子数について

　特異値分解の時と同様に、Matrix Factorizationでもユーザ・アイテム因子行列の因子数は可変です。ただし、特異値分解と異なり、Matrix Factorizationには因子数の上限がありません。特異値分解は解析的な手法であるため分解元の行列の形式に依存していましたが、Matrix Factorizationでは評価値行列の形式にとらわれず任意に設定した因子数でユーザ・アイテム因子行列を学習させることができます（ただし、因子が増えればその分計算リソースも消費することになります）。¹

3. Matrix Factorizationの定式化

　一般的なMatrix Factorizationの定式化について整理します。

　これまでの説明を踏まえると、以下のような式を想定できると思います。

\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
\argmin_{p, q} \sum_{u \in U, i \in I} \frac{1}{2}(r_{ui} - p_{u}q_{i}^T)^2

$\argmin_{p, q}$は、「次の式を最小化する$p$、$q$を求める」という意味です。また、$U$は全ユーザ、$I$は全アイテムを表しています。つまり、すべてのユーザ・アイテムの組み合わせの予測誤差が最小になるようにベクトル$p$、$q$を求めるということになります（実際の計算では「”評価値が観測された”すべてのユーザ・アイテムの組み合わせ」になります）。

　ここで要注意なのが、Matrix Factorizationの学習では、観測された（欠損していない）データのみでパラメータを更新していくため、多く評価しているユーザや多く評価されているアイテムに対して過学習となる可能性が高いということです。したがって、過学習抑制のための正則化項が追加された以下の式が、ベーシックなMatrix Factorizationの式となっています。²

\argmin_{p, q} \sum_{u \in U, i \in I} \{ (r_{ui} - p_{u}q_{i}^T)^2 +\lambda(||p_u||^2_2 + ||q_i||^2_2)\}

$||p_u||^2_2$および$||p_u||^2_2$はL2ノルムの2乗で、例えば$||p_u||^2_2 = \left( \sqrt{p_{u1}^2 + p_{u2}^2 + ...} \ \right)^2$となります。
　各パラメータのL2ノルムの2乗を追加しおり、L2正則化が導入されていることが分かると思います。この式に基づき、SGDやALS（Alternative Least Square）³という手法によってパラメータ$p$、$q$を最適化します。ちなみに、何の断りもなく$\frac{1}{2}$を除外しましたが、最小化を目的とした定式化では定数倍の有無は計算結果に影響しないため、一般的な表記に準拠して除外しています（ただし、アルゴリズムとしては$\frac{1}{2}$倍した目的関数を想定した更新式になっていることも多いです）。

4. まとめ

　以上、計3セクションを通し、Matrix Factorizationの考え方を整理しました。欠損値を補完する必要がなく、さらに予測タスクを解くこともできるという点で、特異値分解よりも汎用的に利用できる行列分解であることが分かったと思います。
　さて、これまでの説明では行列分解を使ったレコメンドを「評価値予測」と「アイテムの序列付け」に限定していましたが、実はこれ以外のレコメンドにも活用できます。例えば、アイテム間の類似度を元にレコメンドする「アイテム間型」のレコメンド手法や、「近傍探索」を利用したレコメンドなどが代表的です。次回は、特異値分解・Matrix FactorizationをPythonで実装し、オープンデータを用いて模擬的にレコメンドを実施してみたいと思います。

参考文献

風間正弘・飯塚洸二郎・松村裕也 (2022) 『推薦システム実践入門―――仕事で使える導入ガイド』,オライリー・ジャパン
Yehuda Koren, Robert Bell and Chris Volinsky(2009), Matrix Factorization Techniques for Recommender System, https://datajobs.com/data-science-repo/Recommender-Systems-[Netflix].pdf

1. 逆に解析的な手法ではないがゆえに、（少なくとも筆者の理解では）「因子が重要度順に並ぶ」ことは保証されません。 ↩

2. 元論文では、$\argmin_{p, q} \sum_{u, i \in K} (r_{ui} - q_{i}^Tp_{u})^2 +\lambda(||q_i||^2 + ||p_u||^2)$ という表記になっていますが、計算していることは同じです。 ↩

3. 「片方のパラメータを定数と見なし、もう片方のパラメータを最小二乗法で最適化する」という作業を交互に繰り返す手法です。今後、Matrix Factorizationの学習アルゴリズムを数式から読み解く記事を執筆したいと考えていますので、その際にALSについても整理する予定です。 ↩