【R】スパースな行列をクラスタリングする

仕事で「超スパースな行列」を扱う必要が出てきまして。生のままで扱うとメモリが足りなくなってしまい、さてどうしたものかと思い悩んでおりました。SASでも扱えないことはないんだけれど、計算にアホほど時間がかかるので、Rで疎行列として扱ったらどうだろうか、と思い至ったわけです。
データセットのイメージとしては、「1,000アイテムを扱っているネットショップの、10,000人の購買者データ」です。購買者IDと買ったアイテムの番号の2列になっていて、各購買者は複数の行で構成されています。それを単純に「アイテムを購入したorしていない」という行列に展開すると、1,000×10,000、つまり一千万の要素を持つ行列になってしまいます。でも、各購買者は平均で3～4アイテムしか購入していないので、データが入る要素はたかだか3万～4万です。だとしたら、疎行列（Sparse Matrix）として扱った方が良い。

というわけでやってみました、というのが今日の話。

なにはなくともデータセット
分析する
- k平均法とspherical-k平均法
- 主成分分析
もうちょっとやりたいこと

なにはなくともデータセット

分析するにはデータセットが必要でございまして。上記したようなイメージのデータセットを作ります。データセットの条件は以下の通りです。

アイテム数は1,000、購買者は10,000人
各購買者を4グループに分け、それぞれがアイテム番号1～250,251～500,501～750.751～1000のアイテムを購入する
各購入者が割り当てられるグループは一様乱数に従う
各購入者が買うアイテムの個数は $rate=0.3$ の指数分布に従う
各購買者が購入するアイテム番号は一様乱数に従う
購入者内でアイテム番号が重複した場合は重複排除する

このような条件に従ったデータセットを作ります。

初期設定

# 初期設定（再現性確保のためseedを設定）
maxnum <- 10000
maxcat <- 1000
grpnum <- 4
maxdev <- maxcat / grpnum 
set.seed(4756);

# 各購入者の購入アイテム数を設定（rate=0.3の指数分布に従う）
cnumf <- function(){
  set.seed(4756)
  round(rexp(maxnum, rate=0.3)+0.5)
}
cnum <- cnumf()

# 各購入者のグループ（1～4）を一様乱数で割り当てる
gnumf <- function(){
  set.seed(4756)
  round(runif(maxnum, min=0.5, max=grpnum + 0.49999))
}
gnum <- gnumf()

イメージ通りのデータになっているか確認します。

hist(cnum)

f:id:wanko_sato:20170326165927p:plain

hist(gnum)

f:id:wanko_sato:20170326165931p:plain

アイテム数は指数分布になっていますし、グループは一様乱数になっています。想定通りなので、次に進みます。

データ本体

結構な量のループをぶん回すので、今回はforeachを使ってみました。参考はこちら。

qiita.com

library(foreach)

# 各購入者毎にseedを設定する
seedsVecF <- function(){
  set.seed(4756)
  round(runif(maxnum)*1000000)
}
seedsVec <- seedsVecF()

# 一様乱数でアイテム番号を割り当てる関数を作る
ccatf <- function(i){
  set.seed(seedsVec[i])
  uni <- runif(cnum[i],
               min = maxdev*(gnum[i]-1),
               max = maxdev*gnum[i]-1)
  uni <- unique(round(uni))
  return(uni)
}

# foreachをぶん回す
ccat <- foreach(i=1:maxnum) %do% ccatf(i)

# 各購入者のアイテム数を計算しなおす（重複排除しているため）
clistC <- foreach(i=1:maxnum) %do% rep(i,length(ccat[[i]]))

# foreachはリスト型でかえってくるのでデータフレームに変換
# 疎行列のNAでない要素に入れるためyn=1を設定
cdfC <- data.frame(id = unlist(clistC),
                  cat = unlist(ccat),
                  yn = 1)

で、データフレームができました。念のため行数を確認します。

> nrow(cdfC)
[1] 38495

行数は38,495です。つまり、1,000×10,000の行列のうち、要素が入るのはおよそ四万程度。密度4%程度の疎行列になります。

疎行列に変換する

次に、このデータフレームを疎行列に変換します。Rにおける疎行列の扱いは次の記事を参考にしました。

d.hatena.ne.jp

Matrixパッケージがあれば疎行列を扱えますので、あらかじめMatrixパッケージをインストールしておいてください。

library(Matrix)

# sparseMtrixにするときに引数i,jに0が入らないようにする
cdfC$cat <- cdfC$cat +1

ctblS <- sparseMatrix(i = cdfC$id,
                      j = cdfC$cat,
                      x = cdfC$yn)

これで先ほどのデータフレームを疎行列に変換できます。iは行、jは列、xは各要素に入る値です。今回は購入の有無だけを問題にしていますので、非ゼロ要素はすべて1になります。たとえば非ゼロ要素に頻度を入れたい場合は、頻度データを別途引っ張ってくる必要があるでしょう。そこはやりやすいようにデータを作ってください。
なお、ここで作った疎行列オブジェクトのサイズはobject.size()関数で調べることができて、

> object.size(ctblS)
467368 bytes

ざっくり0.46MBのサイズとなりました。そのまま展開したときの行列の要素数から考えると、ものすごく小さく圧縮できていることがわかります。

さて、この疎行列、なんとか可視化できないだろうか？と調べてみたところ、

abrahamcow.hatenablog.com

こんな記事を見つけました。おお、これは面白いかもしれない。さっそくやってみよー！

image(ctblS)

ﾜｸﾜｸﾃｶﾃｶ

ｷﾀｰ!!!!

f:id:wanko_sato:20170326174439p:plain

ﾅﾝｼﾞｬｺﾘｬｰ!!!!

まぁ、一様乱数なんで当然なんですけどね・・・少なくともこういう疎行列でも可視化する方法があるのがわかっただけでも収穫でした。

分析する

というわけで、データセットができたのでこれを分析にかけます。一番の目的はクラスタリングです。今回はk平均法を使います。

k平均法とspherical-k平均法

library(cluster)
library(skmeans)

# k平均法とspherical-k平均法を比較
km <- kmeans(ctblS, centers=4)
skm <- skmeans(ctblS, 4, method='pclust', control=list(verbose=TRUE))

突然skmeansなるパッケージが出てきましたが、これはなにかというと、

stats.stackexchange.com

というものらしいです。疎行列を扱えるので便利かなー、と思って使ってみています。ちなみに、通常のkmeansも疎行列に対して実行することができました。なので、性能比較をしてみます。最初に割り当てたグループと、クラスタリングの結果分けられたグループがどれだけ一致しているか、でひかくします。

まずは通常のk平均法から。

> table(gnum,km$cluster)
    
gnum    1    2    3    4
   1 1482    0 1042    0
   2 2448    0    0    0
   3 2536    0    0    0
   4  754 1447    0  291

???
元のグループが行、k平均法によるクラスタリングの結果が列に相当します。1列目がなんだかおかしなことになっています。
では、次にspherical-k平均法の結果。

> table(gnum,skm$cluster)
    
gnum    1    2    3    4
   1    0    0    0 2524
   2    0    0 2448    0
   3    0 2536    0    0
   4 2492    0    0    0

百発百中。
こんだけきれいだと逆に気持ち悪いわ。

でも、なんでk平均法でうまくクラスタリングできなかったんだろう？試しにアイテム購入数とクラスタリング結果をテーブルにしてみました。

> table(cnum,km$cluster)
    
cnum    1    2    3    4
  1  2630    9    0    0
  2  1512  355    0    0
  3  1012  351   42    1
  4   550  276  226    4
  5   387  178  195   11
  6   297  126  135   28
  7   205   73  121   35
  8   144   34   83   36
  9   109   29   59   32
  10   91   10   45   35
  11   68    4   27   26
  12   62    2   26   20
  13   34    0   28   15
  14   27    0   14    8
  15   23    0   11    5
  16   14    0   13    6
  17   14    0    3    6
  18   12    0    5    9
  19   10    0    5    6
  20    5    0    0    2
  21    4    0    0    2
  22    4    0    1    0
  23    2    0    0    1
  24    2    0    1    1
  26    1    0    0    0
  27    0    0    1    1
  29    0    0    0    1
  31    1    0    1    0

k平均法でクラスター1（がちゃがちゃになっていたクラスター）に割り当てられているものが、主に購入アイテム数の少ない購入者であろうことが推定されます。もちろんそれがすべてではないんでしょうけれど、アイテム数が少ないとうまくクラスタリングできないのじゃないかな、と思います。つまるところ、疎行列のクラスタリングではk平均法の性能が劣る、ということがもしかしたらいえるのかもしれません。
※詳しいことはわからないので、教えていただけると嬉しいです。

主成分分析

上の方でimage()関数で無理やりデータを可視化したわけですが、他にうまい可視化の方法がないかな、と思って考えておりましたところ、「そういえば主成分分析があるじゃないか」と思い至りました。思い至ったのでさっそくやってみます。ちなみにprincomp()は疎行列を引数にとることができるようです。

pcaTest <- princomp(ctblS)$scores[,1:2]
plot(pcaTest)

はい、その結果はと言いますと。

f:id:wanko_sato:20170326190857p:plain

当たり前すぎて生きるのがつらい。
そりゃそうだよ、分けた4グループに割り当てたアイテム番号がばっきりわかれてるんだからさ。
だったらはじめから主成分分析しなよ、という話にもなりかねないのでとりあえずこの話はここまで。

もうちょっとやりたいこと

今回作ったデータを眺めていて、改めて思いついたことがありました。

> head(ccat)
[[1]]
 [1] 205 123 245 218  38 176  68 200  55 165 186   8

[[2]]
[1] 487 477 432

[[3]]
[1] 394 326

[[4]]
[1] 215 207 205   1

[[5]]
 [1] 699 724 711 694 523 623 606 747 502 700 664 653 611 659 732 587 706 727 591 573
[21] 576

[[6]]
[1] 270

最初の方で作った各購入者の購入アイテムリストです。list型のオブジェクトに、各購入者のアイテム番号がベクトルで格納されている形になっています。

これって、各アイテム番号をひとつの単語とみなすと、文書集合とみることができて、ということは自然言語処理のスキームをあてはめることができるんじゃね？

と、しょうもないことを思いついてしまったのです。つまり、このデータを半角スペースで区切られた単語集合の集合、つまり文書集合とみなして、形態素解析を行ってDocument-Term-Matrixにして、トピックモデルを当てはめる、とかできるんじゃなかろうか、と。それはそれで、k平均法を当てはめるのとは違うアプローチで面白いんじゃなかろうか、と。

今回は時間の都合でやってませんけれども。