2014-08-23

線形代数の用語と意味まとめ（主に自分用）

恥ずかしながら，線形代数周りの用語って似たようなものが多くて，すぐにアレがどれだっけと混同してしまいがちになります．線形代数の手計算とかがんばってたのなってもう10年とか昔の話だし，チートシート的にまとめなおしておこうと思いました．内容的には，主に統計や機械学習で使うような内容が中心になっています．

概要

統計・機械学習で使う線形代数は，基本的には以下「計算の簡便化」と「データ変換」の2つがメインです．もちろん数学的に突っ込んでいったり，統計・機械学習でも応用的な手法を用いる場合はその限りではないですが，基本的には下の2つが大きいと思います*1．

計算の簡便化

（例えば固有値・固有ベクトルを用いて）行列を対角化することで，行列の乗算を高速に実施する
（LU分解を用いて）扱いやすい形に行列を分解することで，その後の計算を高速にする

データ変換

SVDを行うことでLSIやPCAといったデータ縮約を行うことができる
行列を用いてデータを抽象化することで，回帰係数の算出を統一的な形で記述できる*2

ということで，上のような目的に向けて順に用語と意味を並べました．できるだけ頭から順に読んでいけば理解できるように書き下したつもりです．またあくまで個人的な理解なので，数学的な厳密さについて欠ける箇所や，説明不足な所などあるかもしれません．

個別項目

ノルム

意味: ベクトルの大きさのようなもの
用例: $\vec{a}=(3,4,5)$ のとき $\|\vec{a}\|=\sqrt{3^2+4^2+5^2}$
備考: ↑は2次のノルムで，単に絶対値を加算しただけの1次ノルム，3乗和の3乗根を取った3次ノルムなどいろいろなパターンがある．とはいえ2次ノルム以外は余り使われない*3．

内積

意味: ベクトルの要素同士のかけ算
用例: $\vec{x}=(3,4),\vec{y}=(2,5)$ のときの $\vec{x}\cdot \vec{y}=(\vec{x},\vec{y})=3 \times 2 + 2 \times 5 = 16$
備考: 図形的には $\vec{x}$ から $\vec{y}$ への正射影を意味するので， $\vec{x}$ と $\vec{y}$ のなす角を $\theta$ としたときに， $\|\vec{x}\|\|\vec{y}\|cos\theta$ とも表せる

一次結合

意味: 複数のベクトルを実数倍したものの和
用例: $k\vec{a}+l\vec{b}$

基底

意味: 定義空間 $W$ をベクトルの一次結合で一意に表せる場合に， $\{\vec{a}, \vec{b}\}$ は $W$ の基底という

標準基底

意味: 要素の1つだけが1で，それ以外が0のベクトル群で構成された基底
用例: 定義空間が $R^2$ のときの $\vec{e_{1}}=(1, 0),\vec{e_{2}}=(0, 1)$

正規直交基底

意味: 基底を構成するベクトルがすべてノルム1で，かつすべての要素ベクトルが互いに直交している((要素ベクトル同士の内積がゼロになる)もの
用例: $\vec{e_{1}}=(\frac{2}{\sqrt{5}},\frac{1}{\sqrt{5}}),\vec{e_{2}}=(-\frac{1}{\sqrt{5}},\frac{2}{\sqrt{5}})$ であれば． $|\vec{e_{1}}|=|\vec{e_{2}}|=1$ かつ $\vec{e_{1}} \cdot \vec{e_{2}}= 0$
備考: 標準基底は当然正規直交基底でもある．正規直交基底を用いたベクトル $\vec{a}=a_1\vec{e_{1}}+a_2\vec{e_{2}}, \vec{b}=b_1\vec{e_{1}}+b_2\vec{e_{2}}$ については， $|\vec{a}|=\sqrt{a_1+a_2}$ ， $\vec{a}\cdot \vec{b}=a_1 b_1+a_2 b_2$ と標準基底同様の簡便な計算でノルムや内積が計算できる

線形写像

意味: $n$ 次元ベクトル空間 $R^n$ の1つの要素を定めると，別のベクトル空間 $R^m$ のただ一つの要素が定まる場合，fを線形写像もしくは1次写像という
備考: 特に $f: R^n \rightarrow R^n$ の線形写像については，線形変換もしくは1次変換という．図形的には標準基底で $\vec{OP}=k\vec{e_1}+l\vec{e_2}$ と表される点 $P$ を，基底 $\vec{a},\vec{b}$ で $\vec{OQ}=k\vec{a}+l\vec{b}$ と表される点 $Q$ に移す作業を意味する

単位行列

意味: 対角要素がすべて $1$ で，それ以外がすべて $0$ である行列
用例: $E=\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}$

逆行列

意味: 行列 $A$ に掛けると単位行列 $E$ になる行列
用例: $AA^{-1}=E$

正則

意味: $A$ が逆行列を持つとき $A$ は正則であり，持たない場合は正則ではない

行列式

意味: 行列の面積，体積のようなもの
用例: 2次の場合 $\rm{det} \begin{vmatrix}4&2\\3&7\end{vmatrix}=4\times7-2\times3=22$
備考: $\rm{det}A=0$ のときは， $A$ は逆行列を持たない．なお $|A|$ と書く場合もある*4

線形変換における基底の取替

意味: $A$ による線形変換 $f: X \rightarrow AX$ において，基底を $(\vec{e_1},\vec{e_2}) \rightarrow (\vec{a},\vec{b})=P$ に取替える場合， $A \rightarrow P^{-1}AP$ と表せる

対角行列

意味: 行列の対角要素に0以外の数が，それ以外の要素には0が並ぶ行列のこと
備考: 対角行列は累乗が $\begin{pmatrix}2&0&0\\0&3&0\\0&0&4\end{pmatrix}^n=\begin{pmatrix}2^n&0&0\\0&3^n&0\\0&0&4^n\end{pmatrix}$ と非常に簡単に求められる，といった計算上の扱いやすさがある

固有値と固有ベクトル

意味: $A$ に対して， $A\vec{p}=\lambda\vec{p}$ となるような $\lambda$ を固有値， $\vec{p}$ を固有ベクトルという
備考: 下に示すように固有値と固有ベクトルは，行列を対角化して扱いやすくするための道具と捉えることが可能

固有ベクトルによる対角化

意味: 固有ベクトルを基底とした行列を用いて， $A$ の基底を取替えると，固有値を要素とした対角行列に変換できる
用例: $P=(\vec{p_1},\vec{p_2},\vec{p_3})$ のとき， $P^{-1}AP=\begin{pmatrix}\lambda_1&0&0\\0&\lambda_2&0\\0&0&\lambda_3\end{pmatrix}=L$

スペクトル分解

意味: 対角行列を $E=\begin{pmatrix}1&0\\0&0\end{pmatrix}$ や $E=\begin{pmatrix}0&0\\0&1\end{pmatrix}$ の1次結合で表すこと
用例: 固有値行列 $L$ と固有ベクトルを基底とした行列 $P$ としたとき， $A=PLP^{-1}=\lambda_1 P\begin{pmatrix}1&0&0\\0&0&0\\0&0&0\end{pmatrix}P^{-1}+\lambda_2 P\begin{pmatrix}0&0&0\\0&1&0\\0&0&0\end{pmatrix}P^{-1}+\lambda_3 P\begin{pmatrix}0&0&0\\0&0&0\\0&0&1\end{pmatrix}P^{-1}=\lambda_1 X+\lambda_2 Y+\lambda_3 Z$ と表せる*5
備考: 上の $X, Y, Z$ については， $X^2=X\\Y^2=Y\\Z^2=Z\\XY=YZ=ZX=O\\X+Y=Y+Z=Z+X=E$ という性質がある

転置行列

意味: 行列の列要素と行要素を入れ替えたもの
用例: $B=\begin{pmatrix}2&7&0\\0&3&1\\0&0&4\end{pmatrix}$ のとき， $B^t=\begin{pmatrix}2&0&0\\7&3&0\\0&1&4\end{pmatrix}$

正方行列

意味: 行要素の数と列要素の数が等しい行列．要するに正方形の行列
用例: $B^t=\begin{pmatrix}2&0&0\\7&3&0\\0&1&4\end{pmatrix}$ は $3\times3$ の正方行列である

(実)対称行列

意味: $\begin{pmatrix}2&7&0\\7&3&1\\0&1&4\end{pmatrix}$ のように，行列の上三角成分と下三角成分が対角成分を挟んで対称になっている正方行列
備考:一般的に対称行列は，実数のみを要素として持つ実対称行列を意味することが多い（複素数を要素に持つ場合，対称行列に対応するものとしてエルミート行列がある）．対称行列の場合，転置しても元の行列と等しく $A^t=A$ となる．また対称行列の固有値を $\lambda=(\lambda_1,\lambda_2,\lambda_3)$ ，固有ベクトルを $P=(\vec{p_1},\vec{p_2},\vec{p_3})$ としたとき，異なる固有値に対する固有ベクトルは直交するため， $\vec{p_1}\cdot \vec{p_2}=\vec{p_2}\cdot \vec{p_3}=\vec{p_3}\cdot \vec{p_1}=0$ となる．なお実対称行列の固有値は必ず実数になる*6

直交行列

意味: 行列とその転置行列の積が単位行列になるもの
用例: $CC^t=E$
備考: 対称行列の固有ベクトルを正規化（ノルムが1になるように）した行列を $Q=(\vec{q_1},\vec{q_2},\vec{q_3})$ とすると， $Q^{-1}=Q^t$ が成り立つ．また任意の行列 $A$ について， $AA^t=A^tA$ が成り立つとき， $A$ は直交行列で対角化可能となる．あるベクトル $x$ に直交行列 $A$ をかけてもノルムは変化しない，つまり $\left|Ax\right|=\left|x\right|$ が成り立つ*7

随伴行列

意味: 複素数を成分にとる行列に対して，転置を行った上で複素共役を取ったもの
用例: $\begin{pmatrix}2&2+i\\1&3\\-i&7\end{pmatrix} \rightarrow \begin{pmatrix}2&1&-i\\2+i&3&7\end{pmatrix} \rightarrow \begin{pmatrix}2&1&i\\2-i&3&7\end{pmatrix}$ のように，転置して，複素共役を取った最後の行列は，もとの行列の随伴行列となる

エルミート行列

意味: 複素数要素を持つ正方行列（=複素正方行列）のうち，随伴行列がもとの行列と一致するもの
用例: $B=\begin{pmatrix}2&2+i&3\\2-i&3&i\\3&-i&4\end{pmatrix}$ は $B_{1,2}=2+i$ の対称となる要素が $B_{2,1}=2-i$ となっている
備考: 複素数要素を持たない場合には，単なる対称行列になる．エルミート行列の固有値は必ず実数になる

ユニタリ行列

意味: 複素正方行列のうち，随伴行列が逆行列と等しいもの（=掛け合わせると単位行列になる）
用例: $U=\begin{pmatrix} \frac{1}{\sqrt{2}}&\frac{i}{\sqrt{2}}\\\frac{i}{\sqrt{2}}&\frac{1}{\sqrt{2}} \end{pmatrix}, U^{*}=\begin{pmatrix} \frac{1}{\sqrt{2}}&\frac{-i}{\sqrt{2}}\\\frac{-i}{\sqrt{2}}&\frac{1}{\sqrt{2}} \end{pmatrix}$ のとき， $UU^*=U^*U=I$ より $U^*=U^{-1}$
備考: 直交行列の複素数バージョンと考えればよい．直交行列と同様，ユニタリ行列をかけても対象ベクトルのノルムは変化しない

正定値行列

意味: エルミート行列 or 対称行列 $A$ について，すべての固有値の符号が正であるもの
用例: 固有値 $\lambda=2, 5$ を持つ $\begin{pmatrix}3&\sqrt{2}\\\sqrt{2}&4\end{pmatrix}$ は正定値行列となる
備考: 固有値が0と正の符号で構成されているときは半正定値行列という

三角行列

意味: 対角線の左下，または右上がすべてゼロであるような正方行列のこと．左下がゼロの場合は上三角行列，右上がゼロの場合は下三角行列という
用例: $U=\begin{pmatrix}2&3&5\\0&3&1\\0&0&4\end{pmatrix}$ は上三角行列， $U=\begin{pmatrix}2&0&0\\1&3&0\\4&5&4\end{pmatrix}$ は下三角行列
備考: 三角行列の行列式は，対角成分の積なので非常に簡便に求められる．下のLU分解は，この正方行列を三角行列に分解することで，計算を簡便化するための手法である

LU分解

意味: 正方行列 $A$ を下三角行列 $L$ と上三角行列 $U$ の積に分解すること
用例: $A=\begin{pmatrix}2&4&0\\2&6&4\\3&9&8\end{pmatrix}$ のとき， $A=LU=\begin{pmatrix}1&0&0\\1&1&0\\\frac{3}{2}&\frac{3}{2}&1\end{pmatrix} \begin{pmatrix}2&4&0\\0&2&4\\0&0&2\end{pmatrix}$
備考: 計算しやすくするため $L$ の対角成分を $1$ に揃えるのが慣習となっている．また $L$ も $U$ も対角成分を1にするかわりに，対角行列 $D$ を加えて $A=LDU$ としたLDU分解の形もある．LU分解を行うことで，もとの正方行列を用いた一次方程式が簡単に解けたり，行列式が簡単に求まったりするようになる

コレスキー分解

意味: LU分解の特殊な場合で，正定値エルミート行列 $A$ を下三角行列 $L$ とその随伴行列 $L^*$ に分解すること
用例: $A=LL^*$

QR分解

意味: 行列 $A$ を直交行列 $Q$ と上三角行列 $R$ の積に分解すること
用例: $A=\begin{pmatrix}1&2&3\\-1&0&-3\\0&-2&3\end{pmatrix} = \begin{pmatrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{3}}\\\frac{-1}{\sqrt{2}}&\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{3}}\\0&\frac{-2}{\sqrt{6}}&\frac{1}{\sqrt{3}}\end{pmatrix}\begin{pmatrix}\sqrt{2}&\sqrt{2}&\sqrt{18}\\0&\sqrt{6}&\sqrt{6}\\0&0&\sqrt{3}\end{pmatrix} =QR$
備考: QR分解を行うことによって，固有値を高速に求めることができるようになる

特異値分解(SVD)

意味: $m \times n$ 行列 $M$ を， $m \times r$ の列直交行列 $U$ と $r \times r$ 行列 $\Sigma$ ， $n \times r$ の列直交行列 $V$ に分解すること
用例: $M=U\Sigma V^t$
備考: 列直交とは， $U$ が4行2列のとき $UU^t=\begin{pmatrix}1&0\\0&1\end{pmatrix}$ になることをいう．もちろん $M$ が正方行列であれば， $U,V$ はともに直交行列となる．LSI（潜在意味インデキシング）やPCA（主成分分析）を行う際に使われる

参考

上記のまとめを書くのに使ったのは，主に下記の2冊です．どちらもいわゆる教科書というよりは，線形代数がどう使われるか，何の役に立つのかということを説明している良書です．順番は「意味がわかる線形代数」→「プログラミングのための線形代数」で，前者はいわゆる文系出身で数学忘れました，な人でも理解できるような形で丁寧に書かれている良書です．後者はプログラマ向けで，実装コードや処理の高速化といった面に結構ページを割いています．

まずはこの一冊から意味がわかる線形代数 (BERET SCIENCE)

作者: 石井俊全
出版社/メーカー: ベレ出版
発売日: 2011/06/22
メディア: 単行本
購入: 1人クリック: 5回
この商品を含むブログを見る

プログラミングのための線形代数

作者: 平岡和幸,堀玄
出版社/メーカー: オーム社
発売日: 2004/10
メディア: 単行本
購入: 27人クリック: 278回
この商品を含むブログ (89件) を見る

*1:なので，線形代数の教科書には必ず出てくるジョルダン標準形については，ここでは触れません．あんまり統計・機械学習分野でジョルダン標準形を使う場面をみないし，流れがわかりにくくなるので...

*2:さらには，コレスキー分解やQR分解を用いて最小二乗法を解くことで，回帰分析の係数が定まるといった側面もあります

*3:ブコメをウケて修正しました．

*4:ブコメを受けて修正しました

*5:ブコメを受けて修正しました．

*6:ブコメを受けて追記しました

*7:ブコメを受けて追記しました

2014-08-23

非負値行列因子分解（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))

*1:協調フィルタリングの大枠については，EulerDijkstraさんの記事や神嶌先生の解説論文を参照いただければと思います

*2:それぞれ「千と千尋の神隠し」「となりのトトロ」「おもひでぽろぽろ」「魔女の宅急便」です．

*3:．よくあるAPIのイメージとしては，ユーザIDをパラメタとしてリクエストを投げると，スコアの高い順にソートされた映画IDリストが返ってくる感じです．

*4:t(M)はMの転置行列ですね．

*5:もちろん，自分自身には高いスコアがついてしまいますので，リストから除外しておきましょう．

*6:ちなみにスコアは積和ですので，行と列をひっくり返した「となりのトトロ」と「千と千尋の神隠し」のスコアも同じ20になります．つまり得られた行列は必ず対角行列になります．

*7:ここのところの説明は，細かい意味まで正確かというとあまり自信もないので，イメージとしてこんな感じだ，程度で理解してもらえればと思います．

*8:もちろん，元のMが得られても何の意味もありませんが...

*9:よく似た例が，回帰分析で身長と年齢から体重を予測したときに，値を小さくしていくと体重がマイナスになってしまう，というよくあるパターンと同根の問題だといえます．

*10:全然関係ないですが，同じようにゼロの要素が増えて次元圧縮が行われやすくなるというと，L2正則と比べたときのL1正則を思い浮かべますね．これのグラフィカルな説明はunnonounoさんの記事がわかりやすいです．

*11:そもそも手法が違う時点で，正規化していない限りスコアの絶対値の違いには何の意味もありませんが...

*12:翻訳版が出る前に原書を買ってたのを思い出して，さっき本棚から引っ張りだして読んでました．