#はじめに
こんにちは。@M4_ms1です。この記事は明治大学 Advent Calendar 2018の16日目の記事です。この記事の内容は元々「実験データ解析演習」の最終レポートに出したものだったんですが、今までなんとなく勉強してきた数学をわかりやすい形で応用できたことに感動したので投稿しちゃいます。みんなにもこの感動が伝わるといいな！！（ふんわりした理解のところもあるので（Rの書き方とか特異値分解あたりとか）、間違ってるところがあったらたくさん指摘していただきたいです。）

Rで犬種レコメンダーを作るとは、犬種ごとの情報をまとめた文書群を潜在意味解析(LSA : Latent Semantic Analysis / LSI : Latent Semantic Indexing)して、与える検索質問文に合う犬種を出力させるということです。今回は犬の性格・習性についてまとめた文書を複数用意しました。これらの文書群はジャパンケネルクラブで公開されている「世界の犬」データの、「習性/性格」を利用しました(潜在意味解析は「攻撃的ではない」など否定の文書の意味を解析できないので、そのような文は一文単位で削除したり、同じ意味で書き換えたりした。)。

最終的に、これらの犬種の中から、検索質問文に合う性格・習性の犬種を検索することが目標です。

#潜在意味解析の雰囲気
潜在意味解析では、文書群を語句-文書行列という形で表現します。この行列の要素は（だいたい）ある文書におけるある語句の出現頻度です。この語句-文書行列を特異値分解することで、類似度の高い語句を同じものとみなし、各文書の潜在的な意味をあぶりだします。
例えば、

この本はわかりやすい。
この書籍はわかりやすい。
この本は面白い。
文書1,2は、潜在的には同じ意味です。しかし「わかりやすい本」という検索質問文でこれらが含まれる文書群を検索したとき、「本」と「書籍」は異なる単語なので文書2はヒットしません。そこで、「本」「書籍」という語句は類似語であるとみなすことで、これら2つの文書が同じ意味であるとみなされ、どちらの文書も検索にヒットするようになります。

#語句-文書行列
語句-文書行列$D$は文書群を表現する行列で、行は文書を表し、列は語句を表します。つまり、$n$個の文書、$m$個の語句について語句-文書行列を作成するとき、語句-文書行列$D$は$n\times m$行列になります。各要素は対応する文書における対応する語句の「重み」です。この重みは、その文書でのその語句の出現頻度より計算されます。
例として、各要素は単純に語句の出現頻度であるとし、上の例文で語句-文書行列を作成すると、以下のようになります。

	文書1	文書2	文書3
本	1	0	1
わかりやすい	1	1	0
書籍	0	1	0
面白い	0	0	1

この例では、語句-文書行列の要素を単純に語句の出現頻度としましたが、この出現頻度に重み付けすることで、より検索の精度を上げることができます。各要素の値$D_{ij}$は、局所的重み$l_{ij}$、大域的重み$g_i$、文書正規化係数$n_j$より、以下のように計算されます。
$$D_{ij}=\frac{l_{ij}g_i}{n_{j}}$$
一応各重みにはたくさんの計算方法がありますが、全部書いてるととんでもなく長くなっちゃうので、今回使ったものだけ紹介します。その他のものについては「情報検索アルゴリズム」がわかりやすいです(重み付けに限らず、潜在意味解析全体についてとてもわかりやすく書いてあります)。

語句-文書行列の行ベクトルを語句ベクトル$w_i$、列ベクトルを文書ベクトル$d_j$とします。また、文書$d_j$における語句$w_i$の出現頻度を語句頻度$f_{ij}$、文書の総数を$n$、語句$w_i$を含む文書数を文書頻度$n_i$とします。
##局所的重み
局所的重み$l_{ij}$は語句$w_i$の文書$d_j$における出現頻度を元にして計算されます。出現頻度の高いものほどより重い重みが与えられます。今回は、単純に語句頻度を局所的重みとしました。
$$l_{ij}=f_{ij}$$
##大域的重み
大域的重み$g_i$は、語句$w_i$の出現頻度の偏りを元にして計算されます。全体としての出現頻度が同じでも、多くの文書にまんべんなく出現する語句は特徴が少なく、重要度が低いと考え、逆に限られた文書に集中して出現する文書はその文書の特徴をよく表し、重要度が高いと考えられます。この考えより、後者のような語句により重い重みを与えます。今回は、文書頻度の逆数を採用しました。
$$g_i=\log\frac{n}{n_i}$$
$n_i/n$が小さいほどその語句が出現する文書が少なく、重要であると言えるため、$n_i/n$の逆数をとっています。対数をとっているのは、対数化しないと値の変化が大きすぎるからです。
##文書正規化係数
文書正規化係数$n_j$は、文書の長さによる影響を緩和するための重みです。長い文書に出現する語句は他の語句より出現しやすく、逆も然りです。そのため、各文書にかけられる重みの2乗和が全て1になるように正規化します。これをコサイン正規化といいます¹。
$$n_j=\sqrt{\sum_{i=1}^{m}(l_{ij}g_i)^2}$$

#形態素解析
英語などの文書は単語ごとにスペースで区切られていますが、日本語はそうではないです。当たり前ですが語句-文書行列を作るためには単語ごとの出現頻度を数える必要があり、そのためには単語が単語ごとに区切られている必要があります。このように文書を単語ごとに分ける作業を分かち書きといいます。さらに、上の例では特に断りなく省きましたが、語句-文書行列には「てにをは」などの不要語はカウントしません。日本語の不要語は助詞、助動詞と、「こと」「なる」などのその他の品詞の単語が含まれます。助詞、助動詞以外の不要語を見つけるのを自動で行うのは難しいですが、品詞を特定することはかなり高い精度で自動化できているようです。このように、文書を分かち書きし、さらに品詞を特定する一連の作業を形態素解析といいます。RではRMecabというパッケージが形態素解析に便利です。以下のソースコードは犬種のデータに対して語句-文書行列を作成するまでソースコードです。

library(RMeCab)
category = read.csv("犬種.csv")
res = docMatrix("texts")#語句-文書行列の作成
res.func = res[c(-1,-2,-14,-18,-22,-26,-28,-32,-52,-53,-62,-69,-122,-142,-155),] #手動で機能語だけを抽出
ni = numeric(dim(res.func)[1]) #文書頻度niを全て0に初期化
for(i in 1:dim(res.func)[1]){ #文書頻度を数える
  for(j in 1:dim(res.func)[2]){
    if(res.func[i,j] > 0){
      ni[i] = ni[i]+1
      }
   }
}
gi = log(dim(res.func)[2]/ni) #文書頻度の逆数を設定
nj = numeric(dim(res.func)[2]) #文書正規化係数を全て0に初期化
for(j in 1:dim(res.func)[2]){ #コサイン正規化で文書正規化係数係数を設定
  for(i in 1:dim(res.func)[1]){
    nj[j] = nj[j] + (res.func[i,j]*gi[i])^2
    }
  nj[j]=sqrt(nj[j])
}
res.func.nml = res.func*gi/nj #gi,njを用いて正規化

犬種.csvは、後述するtextsディレクトリに習性・性格が記述してある犬種を、ディレクトリに収められている順番で縦に書き並べたcsvファイルです。3行目のdocMatrix(”texts”)という関数は、ワーキングディレクトリ内にあるtextsというディレクトリの中にあるテキストデータを用いて語句-文書行列を作成する関数です。これらのテキストデータは文書ごとに分かれており、文書の内容だけが書かれている.txtデータになっています。

#特異値分解
$m\times n$の語句-文書行列$D$は、以下のようにエコノミーサイズの特異値分解ができます。
$$D=U\Sigma V\tag{1}$$
ここで、rank($D$)=$r$とすると、$U$は$m\times r$の直交行列、$V$は$r\times n$の直交行列、$\Sigma$は$r\times r$行列です²。$\Sigma$の対角成分は特異値が降順で並んでいて、それ以外の要素は0です。語句-文書行列が特異値分解できることの証明は省きますが、これについても「情報検索アルゴリズム」にわかりやすく書いてあります。
ところで「エコノミーサイズ」って何だよって思った方もいらっしゃると思うので、こちらで解説しました。
特異値分解できて何が嬉しいかというと、$U$と$V$が直交行列なんです。$U$と$V$が直交行列で何が嬉しいかというと、直交行列の全ての行ベクトル（列ベクトル）は正規直交基底を成します。
$D$の列ベクトル$d_ｊ$に注目してみましょう。$A=\Sigma V^T$とおくと、式$(1)$より、$d_j$は以下のように書けます$(1\leq i\leq n,,1\leq j\leq m)$。

\begin{align}
d_j&=Ua_i\\
&=[u_1,\,u_2,\cdots,\,u_n]\,a_i
\end{align}

$U$の列ベクトルは正規直交基底をなすので、$U$は$D$の列ベクトルがなすベクトル空間(=文書ベクトル空間)の正規直交基底いなっていると言えます。同様に考えれば、$V$が語句ベクトル空間の正規直交基底になっていることもわかります。
(この部分の説明はこちらを参考にさせていただきました。ありがとうございます)
さらに、これらの基底は分散が最大になるようにとられていて、$\Sigma$の対角成分である特異値$\sigma_i$は、$U$の$i$番目の列ベクトルや$V$の$i$番目の行ベクトルを基底とみなしたときのその基底の重要度となっています。

##次元圧縮
ここからが本番です。
特異値分解によって得られた行列$U,V,\Sigma$を用いて、語句-文書行列$D$の次元を圧縮します。$U,\Sigma,V$を低ランクで近似することによって$D$を表現することを考えることにしましょう。$U$の列ベクトルは左側にあるほど分散が大きくなり、基底として重要なものであるので、左側の$k$個($k<r$)の列ベクトルを用いて$U$を近似することができます。$U$の左側から$k$個の列ベクトルだけをそのまま抜き出した$m\times k$行列を$U_k$とおきます。同様に考えて、$\Sigma$の左上の$k\times k$だけをそのまま抜き出した$k\times k$行列を$\Sigma_k$、$V$の上から$k$個の行ベクトルだけをそのまま抜き出した$k\times n$行列を$V_k$とおきます。
$$D_k=U_k\Sigma_k V_k\tag{2}$$
$U_k,\Sigma_k,V_k$を式($1$)の$U,\Sigma,V$に代入して得られる$m\times n$行列$D_k$は、語句-文書行列$D$の近似となっています。$U$の列ベクトルの左側$r$列、$V$の行ベクトルの上側$r$行はそれぞれ$D$の列ベクトル、行ベクトルの基底であるため$U_k,V_k$のランクは$k$であり、$\Sigma$は対角行列であるため$\Sigma_k$のランクも$k$です。よって$D_k$のランクは$k$であり、これは式(2)で$D$を低ランクへ近似していることを意味しています。つまり、次元数はそのままで、ランクだけ下げて情報量を減らしているということです。以下のソースコードは、特異値分解を行い、$D_k$を求めるまでのソースコードです。なお、今回は特異値が1以上になるときの添字を$k$としました。

U = svd(res.func.nml)$u #特異値分解する
V = svd(res.func.nml)$v
sigma = svd(res.func.nml)$d
k = 0 #特異値が1以上になる要素数kを初期化
for(i in 1:length(sigma)){ #kを探す
  if(sigma[i] > 1.0){
    k = i
  }
}
Uk = U[,1:k] #U,V,Sigmaの最初のk個だけを抽出
Vk = V[1:k,]
sigma.k = sigma[1:k]
mtrx.sigma.k = diag(sigma.k) #sigma.kは特異値が並んだだけのベクトルなので、対角行列化
Dk = Uk %*% mtrx.sigma.k %*% Vk #DをDkとして再構成

#検索質問文ベクトルの作成
検索される側の文書がベクトル化されているので、検索質問文もベクトル化する必要があります。検索質問文ベクトルについては、単に語句ごとの出現頻度を要素とし、重み付けは行いません（重み付けを行った意味がなくなってしまいます）。また、このとき、検索質問文ベクトルの語句の内容・順番は語句-文書行列と一致している必要があります。以下のソースコードは、検索質問文を得て、検索質問文ベクトルを作るまでのソースコードです。

FAVDOG = readline("( ??? ) な犬が好き！ : ") #好きな犬のタイプを入力
tmp = unlist(RMeCabC(FAVDOG))　#形態素解析する
lst = names(tmp) %in% c("名詞", "動詞", "形容詞") #名詞、動詞、形容詞だけを抽出
FAVDOG = tmp[lst]
q = numeric(dim(res.func)[1]) #検索質問文ベクトルを初期化
q.tmp = rownames(res.func) %in% FAVDOG #語句-文書行列の中にFAVDOGの中身が存在しているかチェック
for(i in 1:dim(res.func)[1]){ #FAVDOGの中身が存在する語句だけ出現頻度をカウント
  if(q.tmp[i]){
    for(j in 1:length(FAVDOG)){
      if(rownames(res.func)[i] == FAVDOG[j]){
        q[i] = q[i]+1
      }
    }
  }
}

なお、検索質問文の処理より前の部分のプログラムは語句-文書行列の処理を行っているだけなので、検索質問文を変えたいときはこのソースコード以降だけを再計算すれば大丈夫です。いろんないぬが気になる欲張りさんにもぴったり。

#類似度の計算
語句-文書ベクトルと検索質問文ベクトルが用意できたので、あとは検索質問文ベクトルとそれぞれの文書ベクトルの類似度を計算し、類似度が上位のものを表示すれば終わりです。この類似度の計算には、コサイン尺度がよく用いられるそうです。文書ベクトル$d_i$と検索質問文ベクトル$q$のコサイン尺度$\cos(d_i,q)$は、以下のように求められます。
$$\cos(d_i,q)=\frac{d_i\cdot q}{||d_j||,,||q||}$$
コサイン尺度は、文書ベクトルと検索質問文ベクトルのなす角のコサインを表しています。コサイン尺度は$0\leq\cos(d_i,q)\leq1$の値をとり、1に近いほど類似性が高いと言えます。以下のソースコードは、コサイン尺度を計算し、類似度が高い3種の犬種を表示するソースコードである。

cos = numeric(dim(res.func)[2]) #コサイン尺度を設定
zettaiti = function(x){ #普通にベクトルのユークリッドノルムを計算するだけの関数が見つからなかったので自作
  y = 0
  for(i in 1:length(x)){
    y = y + x[i]^2
  }
  return(sqrt(y))
}
for(i in 1:dim(res.func)[2]){ #コサイン尺度を計算
  cos[i] = (q %*% Dk[,i]) / (zettaiti(q) * zettaiti(Dk[,i]))
}
cos.tmp = cos * (-1) #order()は昇順で順位付けするので、順番を逆に変えるために符号を反転
odr = order(cos.tmp) #コサイン尺度を順位付け
for(i in 1:3){ #犬種を表示
  print(category[odr[i],1])
}

これで完成です。犬種.csvとtextsディレクトリにあたるデータを用意して、以上のソースコードを出た順に実行すればできると思います。

#実行例
「スタミナがあって強くて優しい」犬が好き！
↓
1.ウエルッシュコーギー・ペンブローク
2.ヨークシャー・テリア
3.ダックスフンド

やっぱコーギーが一番かわいいわあ

これのどこにコサインが関係しているのか知ってる人いたら教えてください笑 ↩
本来の特異値分解の定義では$U$を$m\times m$行列、$V$を$n\times n$行列、$\Sigma$を$m\times n$行列とします。しかしRのsvd()がデフォルトでエコノミーサイズの方法を使っており、特にこだわりもなかったのでその通り特異値分解しました。しかし、潜在意味解析を行う分にはどっちでも大丈夫です。エコノミーサイズの分解と本来の定義での分解については、ここがわかりやすいです。 ↩