"Introduction to Information Retrieval" の１８章に従って SVD(Singular Value Decomposition) と LSI (Latent Semantic Indexing) を試した。
といっても、原理は線形代数で固有値とか計算したことある人にはシンプルな話だし、実装は世にいくつもあるからそれを使えばいい(実装すればわかる、っていうもんじゃあないし)。
むしろわかりたいのは「手触り」なので、R を使って対話的、チュートリアル的に触っていこう。
なお、この記事には間違いその他あるかもしれないが、今日はエイプリルフールなので特に問題ないというのはとてもありがたい。
題材は Example 18.4 のサンプル行列。
| d_1 | d_2 | d_3 | d_4 | d_5 | d_6 | |
| ship | 1 | 0 | 1 | 0 | 0 | 0 |
| boat | 0 | 1 | 0 | 0 | 0 | 0 |
| ocean | 1 | 1 | 0 | 0 | 0 | 0 |
| voyage | 1 | 0 | 0 | 1 | 1 | 0 |
| trip | 0 | 0 | 0 | 1 | 0 | 1 |
まず R の行列にする。
> (C <- matrix(c(1,0,1,1,0, 0,1,1,0,0, 1,0,0,0,0, 0,0,0,1,1, 0,0,0,1,0, 0,0,0,0,1), 5)); [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 1 0 0 0 [2,] 0 1 0 0 0 0 [3,] 1 1 0 0 0 0 [4,] 1 0 0 1 1 0 [5,] 0 0 0 1 0 1
特異値分解(Singular Value Decompositions)
R には SVD の関数があるので、それを呼ぶだけ。
> (decomp <- svd(C)); $d [1] 2.1625010 1.5943824 1.2752903 1.0000000 0.3939153 $u [,1] [,2] [,3] [,4] [,5] [1,] 0.4403475 -0.2961744 -0.5694976 5.773503e-01 -0.24640214 [2,] 0.1293463 -0.3314507 0.5870217 3.330669e-16 -0.72719701 [3,] 0.4755303 -0.5111152 0.3676900 -2.146720e-17 0.61435841 [4,] 0.7030203 0.3505724 -0.1549059 -5.773503e-01 -0.15978815 [5,] 0.2626728 0.6467468 0.4145917 5.773503e-01 0.08661399 $v [,1] [,2] [,3] [,4] [,5] [1,] 0.7486230 -0.2864540 -0.2797116 -3.808531e-16 0.5284591 [2,] 0.2797116 -0.5284591 0.7486230 2.859041e-16 -0.2864540 [3,] 0.2036288 -0.1857612 -0.4465631 5.773503e-01 -0.6255207 [4,] 0.4465631 0.6255207 0.2036288 2.285092e-16 -0.1857612 [5,] 0.3250960 0.2198798 -0.1214672 -5.773503e-01 -0.4056409 [6,] 0.1214672 0.4056409 0.3250960 5.773503e-01 0.2198798
これは元の行列を以下のように分解しているはず。
Σは特異値を大きい順に並べた対角行列であり、svd() 関数の返値の $d がその特異値。IIR の Σにちゃんと一致している。
U と V は $u と $v に対応するが、IIR に書かれている値と比べると、いくつかの列で符号が逆転している。単位行列化する際に±の自由度があるので、本質的にはちゃんと同じものである。
さて、これが本当に元の行列の分解になっているかどうかを確かめるには以下の式を試せばいい。
> decomp$u %*% diag(decomp$d) %*% t(decomp$v); [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.000000e+00 -3.705261e-17 1.000000e+00 1.216305e-16 -1.856527e-16 2.205487e-16 [2,] -7.918742e-17 1.000000e+00 9.400033e-17 2.842338e-16 5.648693e-17 3.110169e-16 [3,] 1.000000e+00 1.000000e+00 2.290919e-16 -1.312088e-16 -1.432367e-16 2.590566e-17 [4,] 1.000000e+00 1.354728e-16 -1.734520e-16 1.000000e+00 1.000000e+00 -7.392988e-17 [5,] -9.275519e-17 -1.037658e-16 3.993150e-16 1.000000e+00 4.625562e-17 1.000000e+00
計算誤差の範囲で、ちゃんと元の行列 C と一致していた。
それから、U と V は正規直交基底になっているはず。
> t(decomp$u) %*% decomp$u [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 -2.768781e-17 1.719612e-16 1.113476e-16 -1.807450e-16 [2,] -2.768781e-17 1.000000e+00 -3.976311e-17 1.941535e-16 -1.607262e-16 [3,] 1.719612e-16 -3.976311e-17 1.000000e+00 3.117081e-19 -3.032513e-16 [4,] 1.113476e-16 1.941535e-16 3.117081e-19 1.000000e+00 -2.377215e-16 [5,] -1.807450e-16 -1.607262e-16 -3.032513e-16 -2.377215e-16 1.000000e+00 > t(decomp$v) %*% decomp$v [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 7.930261e-17 -2.528766e-16 6.198926e-17 1.069633e-17 [2,] 7.930261e-17 1.000000e+00 -1.045577e-16 -7.665309e-17 9.616873e-17 [3,] -2.528766e-16 -1.045577e-16 1.000000e+00 9.120851e-17 1.929880e-16 [4,] 6.198926e-17 -7.665309e-17 9.120851e-17 1.000000e+00 -1.921748e-17 [5,] 1.069633e-17 9.616873e-17 1.929880e-16 -1.921748e-17 1.000000e+00
確かに。
ちなみに U の行は CC^T の固有ベクトル。つまり PRML(パターン認識と機械学習) の 12章、主成分分析の固有ベクトルと同じもの。
PRML の図 12.3 に図示されているのがその固有ベクトルであり、ってことは、こいつらは正規直交基底なのか〜と思って眺めるとまた楽しめるかもしれない。
低ランク近似(Low-Rank Approximations)
この分解において、固有値を小さい方から捨てていくことで低ランク近似を行う。
U と V が小さくなって扱いやすくなる、というだけでなく、synomymy(類義性) や polysemy(多義性) の問題を解決できる、かもしれない、というメリットがある、らしい。エイプリルフールなので詳細は省略。
次の式でランク２に落とした C_2 を求められる。
> (C_2 <- decomp$u[,1:2] %*% diag(decomp$d[1:2]) %*% t(decomp$v[,1:2])); [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.8481456 0.5159023 0.28162515 0.12986018 0.20574267 -0.07588249 [2,] 0.3607778 0.3575076 0.15512454 -0.20565325 -0.02526436 -0.18038889 [3,] 1.0032701 0.7182854 0.36077778 -0.05052871 0.15512454 -0.20565325 [4,] 0.9780058 0.1298602 0.20574267 1.02853450 0.61713858 0.41139591 [5,] 0.1298602 -0.3860421 -0.07588249 0.89867432 0.41139591 0.48727840
低ランク「近似」って言うくらいだから、なんか元の行列と多少は似てるんかなあなんて漠然とした印象を持っていたんだけど、この C_2 は元の C とは全く似ても似つかない。
ま、そりゃランク下がってるし、大量の情報捨ててるんだから、似てる方がおかしいんだけど。
> qr(C)$rank [1] 5 > qr(C_2)$rank [1] 2
うん、確かに下がってる。
って、あれれ？ よく見ると IIR に載ってる C_2 と全然違う。
そもそも IIR で C_2 だといっている行列は、下３行が全部ゼロになってて、明らかにおかしい。
これってもしかして、 Σ_k * V_k^T じゃあないの？
> diag(decomp$d[1:2]) %*% t(decomp$v[,1:2]) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.6188981 0.6048766 0.4403475 0.9656932 0.7030203 0.2626728 [2,] -0.4567172 -0.8425659 -0.2961744 0.9973192 0.3505724 0.6467468
ほら、やっぱりー。
C_k と C がどのくらい「似ているか」という話が IIR の Theorem 18.4 。
C_k と元の行列 C とのフロベニウス・ノルム が k+1 番目に大きい固有値の平方根(=Σの k+1 番目の値)になり、それは rank k の近似行列のなかでもっとも最小となる。
よし確かめよう。
> sqrt(sum((C-C_2)^2)) [1] 1.667793 > decomp$d[3] [1] 1.275290
おかしいな……計算誤差のレベルじゃあないし。
あれこれ試したり確認したが、まだちょっと原因わからず。
Th18.4 の証明追いかけてみるかなあ。
【追記】
やっぱり IIR の Th18.4 間違ってるっぽいなあ。
言語と計算 (4) 確率的言語モデル には以下のように書いてある(記号は IIR にあわせた)。
例で計算してみる。
> sqrt(sum((C-C_2)^2)) [1] 1.667793 > sqrt(sum(decomp$d[3:5]^2)) [1] 1.667793
ＯＫ。
無料で読める PDF の方を見ているので、書籍だとちゃんと直っているのかな？
潜在意味インデックス(Latent Semantic Indexing)
IIR のサンプルでは、図を書いて LSI 終わり。えー。
そこ、もうちょっとあれこれ遊ぶところでしょう！
クエリー [ocean, trip] に対して、類似ワード検索してみる。
最初の単語-ドキュメント行列の定義を見ると、それぞれのインデックスは３と５とわかるから、
> q <- c(0, 0, 1, 0, 1); > (q_k <- diag(1/decomp$d[1:2]) %*% t(decomp$u[,1:2]) %*% q); [,1] [1,] 0.34136544 [2,] 0.08506838
これと、V_k の各列との余弦類似度を取れば関連語を考慮した検索となる(ことが期待できる)。
V_k の代わりに復元式に放り込めば、Bayesian Sets みたいなこともできたりするんかな！？
> decomp$u[,1:2] %*% diag(decomp$d[1:2]) %*% q_k [,1] [1,] 0.28489529 [2,] 0.05052871 [3,] 0.28171457 [4,] 0.56652045 [5,] 0.28162515
[ocean, trip] に関連のある単語は、最大のインデックス４にあたる "voyage" だ、という結果。
このデータサイズではなんともいえないけど、一応なんとなくそれっぽく。
