非負値行列因子分解（NMF）によるレコメンドのちょっとした例

最近線形代数についていろいろ読みなおしたりしてるのですが（線形代数チートシートを前の記事でまとめてあります），その一環でレコメンドアルゴリズムについていくつか試してみたので，それを解説します．順序としては，基本の協調フィルタリング（ユーザベースド，アイテムベースド）→特異値分解（SVD）→非負値行列因子分解（NMF）になります．

基本的な考え方

ここで取り扱うのは，すべて以下のようなユーザ×商品のマトリックスをベースとしたレコメンドになります*1．ここでは映画レンタルサービスを例にして考えます．6人のユーザが，4つの映画*2のうちレンタル視聴したものについては，1-5点の5段階評価を行いました．0になっているものは「みていない」ということになります．

まずはざっと評価の状況をみると，「千と千尋の神隠し」が最もよく視聴されていて，6人中4人がみています．次にみられているのは「となりのトトロ」と「魔女の宅急便」が同数で3人．またこの2作品は3人中2人がどちらも視聴しており，評価も高めなのがみて取れます．そして「おもひでぽろぽろ」は1人しか視聴しておらず，また評価も低いことがわかります．

> # 6ユーザが4アイテムを評価
> user_item = c(5,0,0,0,
+               3,4,0,0,
+               2,0,1,0,
+               2,4,0,3,
+               0,5,0,4,
+               0,0,0,5)
> items = c('sen', 'totoro', 'omohide', 'majyo')
> users = c('user1', 'user2', 'user3', 'user4', 'user5', 'user6')
> # ユーザ-アイテムマトリックス
> M  = t(matrix(user_item, 4, 6))
> colnames(M) = items
> rownames(M) = users
> print(round(M, 2))
      sen totoro omohide majyo
user1   5      0       0     0
user2   3      4       0     0
user3   2      0       1     0
user4   2      4       0     3
user5   0      5       0     4
user6   0      0       0     5

サービス提供側としては，ユーザに対して未視聴のおすすめコンテンツを提示して，もっと視聴してもらいたいと考えます．効率よくお勧めするために，ユーザ毎に未視聴のおすすめコンテンツリストが返ってくるようにしたい訳ですが*3，そのためには全コンテンツについておすすめスコアを定義する必要が出ます．このおすすめスコアをできるだけ精度よく推定したい，というのが解決すべき問題になります．

協調フィルタリング

アイテムベースド

最初から例外でアレですが，これだけは，「ユーザ→アイテムリスト」ではなく，「アイテム→アイテムリスト」の形になります．Amazonの「この商品を買った人はこんな商品も買っています」が典型例ですね．これを行列演算で表すと，以下のようになります．

> # アイテムベースのレコメンドスコア算出
> # t(M)とMの内積は，以下のループ処理をまとめて行っているのと同等
> # R_item = matrix(0.0, 4, 4)
> # for (i in 1:4) {
> #   I = M[, i]
> #   R_item[i,] = (I %*% M)
> # }
> R_item = t(M) %*% M
> print(round(R_item, 2))
        sen totoro omohide majyo
sen      42     20       2     6
totoro   20     57       0    32
omohide   2      0       1     0
majyo     6     32       0    50

ここでは t(M)とMの内積 *4を取ることで，アイテム同士の評価値の積和を得ることができます．積和が大きいほど評価値が大きい=おすすめ度合いが高い，と考えられるので，これを降順に並び替えれば，あるアイテムに対するおすすめアイテムリストを得ることができます*5．具体的な計算としては，「千と千尋の神隠し」と「となりのトトロ」の評価値の積和は以下のように計算できます．

$score_{sen, totoro}=\sum_{i=1}^{6} s_{sen, i} s_{totoro, i}=5\times0+3\times4+2\times0+2\times4+0\times5+0\times0=20$

これがアイテム行列の $R\_item[1, 2]$ になるわけです*6．

ユーザーベースド

さて，続いてユーザのほうにいきますが，これも基本的にはアイテムと同じで，今度は順番をMとt(M)の内積を取る形になります．ただ，アイテムベースと違うのは，単に内積を取るだけだと，ユーザに対してユーザをお勧めする形になってしまって，肝心のアイテム推薦ができません．そこで算出したユーザ同士の行列とMの内積を取ることで，類似ユーザスコアに応じたアイテムの評価値の積和を得ることができ，これがユーザID→アイテムリストを得ることができました．

> # ユーザベースのレコメンドスコア算出
> # こちらもアイテムベースと同じ
> # ただし最終的にアイテムをレコメンドするためには，
> # Mとt(M)の内積で類似ユーザを算出した後，
> # 類似ユーザの好むアイテムを得るために，再度Mとの内積を取る必要がある
> R_user_tmp  = M %*% t(M)
> R_user      = M %*% t(M) %*% M
> print(round(R_user_tmp, 2))
      user1 user2 user3 user4 user5 user6
user1    25    15    10    10     0     0
user2    15    25     6    22    20     0
user3    10     6     5     4     0     0
user4    10    22     4    29    32    15
user5     0    20     0    32    41    20
user6     0     0     0    15    20    25
> print(round(R_user, 2))
      sen totoro omohide majyo
user1 210    100      10    30
user2 206    288       6   146
user3  86     40       5    12
user4 182    364       4   290
user5 124    413       0   360
user6  30    160       0   250

この場合，ユーザ1におすすめする動画は，視聴済みの「千と千尋の神隠し」をのぞくと，「となりのトトロ（100）」「魔女の宅急便（30）」「おもひでぽろぽろ（10）」の順番になります．

特異値分解（SVD）

さて，上のシンプルな協調フィルタリングには大きな問題があります．それは．疎な行列に対してうまくレコメンドができないことです．ほとんど視聴した人がいないような映画や，また逆にほとんどサービスを利用していないユーザに対しては，そもそも他のユーザ，アイテムとの視聴の重なりが少ないため，適切なレコメンドを行うことができません．

この問題に対応するための手法の一つが特異値分解（SVD）による行列分解を行うことです．この手法自体の意味については，潜在意味分析（LSI）についての記述ではありますが，あらびきさんの記事がわかりやすく参考になると思います．ここでは一言でざっくりと説明すると，行列のランク数を削減することによって，実質的に似たユーザをクラスタリングして扱うことができるようになります．それによって，ほとんど何のコンテンツもみてないユーザであっても，似たクラスタのユーザの情報をうまく参照することができ，その結果としてアイテムをレコメンドすることができるようになる，という感じかと思います*7．

では，実際にSVDを使ってユーザ→アイテムリストのレコメンドを行ってみると，以下のようになります．ここではランクを2個削減してみました．下から2行目でd_rを使っているところをdを使うように直せば，実際にはランク削減は起こらないので，元の行列と全く同じものが得られます*8．

> # SVD
> res.svd = svd(M)
> u   = res.svd$u
> v   = res.svd$v
> d   = diag(res.svd$d)
> d_r = d
> for (i in 3:4) { # データ圧縮のためランクを落としたdを作る
+   d_r[i, i] = 0
+ }
> R_svd = M %*% v %*% solve(d) %*% d_r %*% t(v)
> # R_svd = M %*% v %*% solve(d) %*% d %*% t(v) #もとのMに一致する
> colnames(R_svd) = items
> print(round(R_svd, 2))
        sen totoro omohide majyo
user1  3.65   1.59    0.82 -1.32
user2  3.79   3.27    0.88  0.62
user3  0.16   0.08    0.04 -0.05
user4  1.94   3.92    0.51  3.06
user5  0.53   4.20    0.21  4.65
user6 -1.32   1.88   -0.24  3.47

ちなみにSVDの場合，uもvも元の行列の直交基底を取っているため，直交行列になっています．SVDというのは，情報を圧縮しているため主成分分析と似たような変数要約的な意味合いを持つ処理です．主成分分析の主成分同士が直交するのと同様，SVDで圧縮された次元も直交します．しかし実際の商品から構成されるユーザ-商品の潜在的な主成分同士が直交するというのは，割と強引な仮定だったりします．

今回の例だと，ユーザと商品から抽出される主成分は，ジャンルと考えることができます．ユーザ視点からみると好みのジャンル，商品視点だと商品の所属ジャンル，という話になります．この場合，「となりのトトロ」と「魔女の宅急便」は「少女の成長ドラマ」というジャンル分けができ，user4，user5がそういったジャンルを好む，とみなすことができます．しかしここで「おもひでぽろぽろ」は，「オトナの人間ドラマ」というジャンルに属知っていると考えられますが，この2ジャンルが直交する（=全く相関がない）というのは，ややしっくりきません．「カーアクション」と「恋愛ドラマ」なら直交といってもいいかもしれませんが，データ特性によっては，直交性の仮定はかなり厳しいものになります．

非負値行列因子分解（NMF）

それでは最後に非負値行列因子分解（NMF: Non-negative Matrix Factorization）を試してみましょう．手法の意味などについては，こちらもあらびきさんの記事がわかりやすくまとまっているので読んでいただければ理屈はわかるかと思います．NMFの特徴は，そのものズバリ分解した行列の要素がすべて正であるという点です．SVDはどうだったかというと，上の例をみていただければわかるように値が負の要素がいくつか存在しています．

この負の要素が存在することの何が問題かというと，実際の評価値は1-5点の正の値しかとらず，また未視聴であっても0点であり，負の値を取り得ないことです．その人へのおすすめスコアといいつつも，実際は0の部分の値を他の人の評価値から推測しているわけなので，それが負の値を取るというのは，論理的におかしなことになります*9．なので，0の値を取らない方が論理的に妥当だといえます．

また，負の要素がないため，積和でスコアを算出する際に，引き算は使えず足し算しか使えません．複数の人の評価の足し合わせでスコアを考えるとして，すべて足し算しか使えない（=凝ったやり方ができない）という制約によって，複雑な解を取りにくくなる方向に圧力が働きます．その結果，引き算が使えるときと比べて，ぐっと結果がシンプルになります．行列演算の文脈でいうと，非ゼロ要素の数が減るという結果として表れます*10．

また直交性についても，NMFでは直交性の仮定が存在しないので，その点においてもSVDより優れていると考えられます．

ということで，これだけは{NMF}パッケージを使って計算します．nmf()の第2引数で圧縮する次元数を指定します．ここではSVDのときと同じ2次元への圧縮を選択しました．またnmf()は初期値の与え方によって毎回異なる結果が得られるので，ここでは適当に1234でseedを固定しています．結果をみるとわかるように，0要素の数が増え，また負の要素がすべてなくなりました．

> # NMF
> library(NMF)
> res.nmf = nmf(M, 2, seed=1234)
> w   = basis(res.nmf)
> h   = coef(res.nmf)
> h_z = rbind(h, rep(0, 4))
> R_nmf = w %*% h
> print(round(R_nmf, 2))
       sen totoro omohide majyo
user1 2.77   1.40    0.83  0.00
user2 4.98   2.52    1.49  0.00
user3 0.55   0.28    0.17  0.00
user4 1.70   3.22    0.51  3.57
user5 0.00   3.58    0.00  5.42
user6 0.00   1.99    0.00  3.01

今回のデータだと，そもそもユーザ-アイテム行列がスパースではないため，眼に見えるほどのおおきな変化はみられず，どの手法においても（スコアの絶対値は違いますが*11），アイテムリストの順番はどれでもそんなに変わらない結果が得られているかと思います．

ちなみにO'Reillyの集合値プログラミングには，Python実装でのNMFのアルゴリズムが載っていますね．中身がどうなっているか気になる方は，そちらを参照しても良いかと思います*12．

集合知プログラミング

作者: Toby Segaran,當山仁健,鴨澤眞夫
出版社/メーカー: オライリージャパン
発売日: 2008/07/25
メディア: 大型本
購入: 91人クリック: 2,220回
この商品を含むブログ (277件) を見る

まとめ

実際のロジックにおける最適化関数

今回の例ですが，あくまでレコメンドアルゴリズムを行列演算の観点から理解するためのサンプルですので，実際にレコメンドエンジンに組み込まれているロジックと厳密には異なっている場合があります．特にMatrix Factorizationの文脈では，SVDもNMFも，CAのHattoriさんの記事にあるように最適化する目的関数は，単なる予測評価値の二乗誤差の最小化ではなく，いくつかのパラメタが追加されています．

上の記事の例だと，アイテム毎，ユーザ毎の平均評価点を考慮した形で二乗誤差を求める形になっています．この最適化問題を解くことによって，やたらと低い点ばかりつけるユーザ，また人気作でみんなが高評価をつける商品を考慮したスコア付けにすることが可能です．そして上の例は，疎な行列の問題にも対応できているのが優れている点です．疎な行列とは何かというと，実際のユーザ-アイテム行列は，ほとんどの要素が未評価（=0）であるということです．

疎な行列への対応

今回の例では，6人×4商品だったので，割と密度の高い行列になっていました．しかし実サービスでのレコメンドを考えると，登録ユーザ数が3万人いて，扱う商品が6000個だったとすると，ほとんどのユーザは，数個程度しか購入しておらず，従ってほぼすべての要素が0で埋め尽くされた，非常に疎な行列になってしまいます．

ここで問題なのは，0というのは評価点0なのではなく，あくまで欠損値だということです．今回の例では，理解しやすくするために，意図的にこのあたりをごまかして説明しました．しかし実際は，評価点0でも実際にみたら4点や5点がつくような商品はいくらでもあるわけです．そうなると，ほとんど0の行列という時点で，二乗誤差の最小化を行うと，この0に大きく引きずられた結果が得られてしまう訳です．

そのため実用的なアルゴリズムの場合，実際に評価された要素だけを対象にして二乗誤差を求めます．もちろん0要素が少ないのであれば，欠損値補完を使ってあげてもいいのかもしれませんが，実際問題として大半が0である場合には，補完も何も...という話なわけです．

追記

ブコメでマイナスのコメントとプラスのコメント両方いただきましたが，個人的にはどちらもその通りだと思っています．理論的な部分はあらびきさんにおんぶに抱っこなのが事実の一方で，このエントリではハンズオンでレコメンドの行列演算を実感できることを目指しています．反応いただけることはありがたいので，今後も反応いただけるような内容のものを出せるよう精進したいです．

最後に，今回のコード一覧はこちらになります．

 # 6ユーザが4アイテムを評価
user_item = c(5,0,0,0,
              3,4,0,0,
              2,0,1,0,
              2,4,0,3,
              0,5,0,4,
              0,0,0,5)
items = c('sen', 'totoro', 'omohide', 'majyo')
users = c('user1', 'user2', 'user3', 'user4', 'user5', 'user6')
# ユーザ-アイテムマトリックス
M  = t(matrix(user_item, 4, 6))
colnames(M) = items
rownames(M) = users
print(round(M, 2))

# アイテムベースのレコメンドスコア算出
# t(M)とMの内積は，以下のループ処理をまとめて行っているのと同等
# R_item = matrix(0.0, 4, 4)
# for (i in 1:4) {
#   I = M[, i]
#   R_item[i,] = (I %*% M)
# }
R_item = t(M) %*% M
print(round(R_item, 2))

# ユーザベースのレコメンドスコア算出
# こちらもアイテムベースと同じ
# ただし最終的にアイテムをレコメンドするためには，
# Mとt(M)の内積で類似ユーザを算出した後，
# 類似ユーザの好むアイテムを得るために，再度Mとの内積を取る必要がある
R_user_tmp  = M %*% t(M)
R_user      = M %*% t(M) %*% M
print(round(R_user_tmp, 2))
print(round(R_user, 2))

# SVD
res.svd = svd(M)
u   = res.svd$u
v   = res.svd$v
d   = diag(res.svd$d)
d_r = d
for (i in 3:4) { # データ圧縮のためランクを落としたdを作る
  d_r[i, i] = 0
}
R_svd = M %*% v %*% solve(d) %*% d_r %*% t(v)
colnames(R_svd) = items
print(round(R_svd, 2))

# NMF
library(NMF)
res.nmf = nmf(M, 2, seed=1234)
w   = basis(res.nmf)
h   = coef(res.nmf)
h_z = rbind(h, rep(0, 4))
R_nmf = w %*% h
print(round(R_nmf, 2))