分からんこと多すぎ

何故か就職できたので技術磨く。一人前の多重線形代数使いを目指しつつ、機械学習、データマイニングをやる

グラフ(ネットワーク状)構造のクラスタリング(グループ分け)基礎:The Markov Cluster Algorithm(MCL)まとめた

MCLのスライド

an introduction to the MCL algorithm

めちゃくちゃ基本らしい。けど微妙に関係ないので、さわりだけ読む。

無意味にグラフ理論の基礎に触れつつ、エッジリストから高速にクラスタを組む話をする。

元論文は有料だったので諦めた。ので、何か微妙に間違っているかもしれない。コードあるけどR言語です。

グラフ構造とは


グラフ構造とは、網目状の構造である。例えば、鉄道とか道路とかネットとか神経回路がそれに相当する。
例えば、道路の例を考えよう。

f:id:rishida:20130606184802p:plain

丸(ノード)が交差点で、線(エッジ)が道路になる。

この図は、「1」の交差点から、2,3,4の交差点にいける事を意味している。

1から5,6,7には行けないし、1から1に行こうと思ったら、1→2→1のように移動しなければならない。

ちなみにこれは無向グラフと呼ばれ、一方通行とかはないものとしている。
一方通行が現れると、有向グラフと呼ばれる。

で、このグラフは行列的に、以下のように表される。

data <- matrix(0,7,7)
data[1,] <- c(0,1,1,1,0,0,0)
data[2,] <- c(1,0,1,1,1,0,0)
data[3,] <- c(1,1,0,1,0,0,0)
data[4,] <- c(1,1,1,0,0,0,0)
data[5,] <- c(0,1,0,0,0,1,1)
data[6,] <- c(0,0,0,0,1,0,1)
data[7,] <- c(0,0,0,0,1,1,0)

> data
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    1    1    1    0    0    0
[2,]    1    0    1    1    1    0    0
[3,]    1    1    0    1    0    0    0
[4,]    1    1    1    0    0    0    0
[5,]    0    1    0    0    0    1    1
[6,]    0    0    0    0    1    0    1
[7,]    0    0    0    0    1    1    0

#冒頭の図
library('igraph')
plot(graph.adjacency(data>0,mode='undirected'),layout=layout.circle)

マルコフ連鎖してみる

とりあえずこの行列をXと呼ぶ。

X_{ij}は、jから出発してiにいけるのかを表している。

行列の各要素の値を、その要素が含まれている列の和で割ると、つまりjを固定してX_{i j}/ \sum_{i} X_{i j}をすると、jから出発してiに行く確率になる。

なお、移動は完全にランダム(等確率)で行うと仮定する。

#列和が1になるように正規化(移動は等確率を成立させる。詳細後述)
data <- t(t(data)*(1/rowSums(data)))

> round(data,2)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.00 0.25 0.33 0.33 0.00  0.0  0.0
[2,] 0.33 0.00 0.33 0.33 0.33  0.0  0.0
[3,] 0.33 0.25 0.00 0.33 0.00  0.0  0.0
[4,] 0.33 0.25 0.33 0.00 0.00  0.0  0.0
[5,] 0.00 0.25 0.00 0.00 0.00  0.5  0.5
[6,] 0.00 0.00 0.00 0.00 0.33  0.0  0.5
[7,] 0.00 0.00 0.00 0.00 0.33  0.5  0.0

行列Xの各列を、それぞれの列和で割ってみた。

これは、移動を等確率にすることを意味する。

例えば1番の交差点からは、2,3,4番の交差点にいける。

(移動を等確率にしたい)=(1→2,1→3,1→4に向かって、それぞれ33%で移動したい)

つまり、3つ道があるなら、3で割れば良い。4つなら4、5つなら5で割れば良い。

したがって、列和で列を割ることにする。

この結果として、X_{ij}には、jからiに移動する確率が入る。

この確率による移動の連鎖が、一次マルコフ連鎖という代物である。

最初に各交差点にいる確率(今回は単位行列を暗に仮定)と、iにいる人が、jに移動する確率は定まっている(マルコフ連鎖

そしてその移動確率は、今いる場所にだけ依存する(一次)

グラフにおけるグループ

さて、X_{ij}は、jからiに移動する確率であるとした。

ここで、行列X同士の積を取る。

X^2とかすると、X^2_{ij} = \sum_k X_{ik} X_{kj}

つまり、X_{ik} X_{kj}は、jからkに行ってkからiに行く確率になるので、結局X^n_{ij}は、jから出発してなんだかんだあってiに到着する確率になる。

今回の場合は消滅したりしないので、jから出発した人はどこかしらにいる。

つまり、X_{ij}のあるjに関して、全てのiについて足し合わせると(行列の列和を取ると)、総和は1になる。

X^2X^3も列和は常に1になる。

> colSums(data%*%data)
[1] 1 1 1 1 1 1 1
> colSums(data%*%data%*%data)
[1] 1 1 1 1 1 1 1

で、図を見ると分かるが(つまりX行列(data)を見ると分かるが)、1,2,3,4がひとつのグループ、5,6,7がもうひとつのグループに見える。

f:id:rishida:20130606184802p:plain:w300

このグループを単純な数値計算で簡単に見つけるのが、MCLである。

MCLは、行列Xを無限回遷移させると(無限回の行列積をする=無限回移動をする)と、グループ内での移動が頻発し、結果として遷移確率が偏ることを想定している。

グループ間にある繋がりは少ないので、グループ間を移動することは、あんまりないと考えるのは自然だろう。

例えば1,2,3,4と5,6,7が、完全に繋がってなければ、無限回移動しても(無限回行列Xの積をとっても)グループ間の移動確率(行列の右上と左下)は0のままである。

#初期化
data <- matrix(0,7,7)
data[1,] <- c(0,1,1,1,0,0,0)
data[2,] <- c(1,0,1,1,1,0,0)
data[3,] <- c(1,1,0,1,0,0,0)
data[4,] <- c(1,1,1,0,0,0,0)
data[5,] <- c(0,1,0,0,0,1,1)
data[6,] <- c(0,0,0,0,1,0,1)
data[7,] <- c(0,0,0,0,1,1,0)

#切り離す
data[2,5]<-0
data[5,2]<-0

#正規化
data <- t(t(data)*(1/rowSums(data)))
tmp <- data

#100回遷移
for(i in 1:100){
  tmp<-tmp%*%data
}

> tmp
     [,1] [,2] [,3] [,4]      [,5]      [,6]      [,7]
[1,] 0.25 0.25 0.25 0.25 0.0000000 0.0000000 0.0000000
[2,] 0.25 0.25 0.25 0.25 0.0000000 0.0000000 0.0000000
[3,] 0.25 0.25 0.25 0.25 0.0000000 0.0000000 0.0000000
[4,] 0.25 0.25 0.25 0.25 0.0000000 0.0000000 0.0000000
[5,] 0.00 0.00 0.00 0.00 0.3333333 0.3333333 0.3333333
[6,] 0.00 0.00 0.00 0.00 0.3333333 0.3333333 0.3333333
[7,] 0.00 0.00 0.00 0.00 0.3333333 0.3333333 0.3333333

だがしかし、元々の行列を無限回遷移させると、

data <- matrix(0,7,7)
data[1,] <- c(0,1,1,1,0,0,0)
data[2,] <- c(1,0,1,1,1,0,0)
data[3,] <- c(1,1,0,1,0,0,0)
data[4,] <- c(1,1,1,0,0,0,0)
data[5,] <- c(0,1,0,0,0,1,1)
data[6,] <- c(0,0,0,0,1,0,1)
data[7,] <- c(0,0,0,0,1,1,0)

#列和で1になるように正規化
data <- t(t(data)*(1/colSums(data)))

tmp <- data
for(i in 1:100){
  tmp<-tmp%*%data
}

> round(data,2)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.00 0.25 0.33 0.33 0.00  0.0  0.0
[2,] 0.33 0.00 0.33 0.33 0.33  0.0  0.0
[3,] 0.33 0.25 0.00 0.33 0.00  0.0  0.0
[4,] 0.33 0.25 0.33 0.00 0.00  0.0  0.0
[5,] 0.00 0.25 0.00 0.00 0.00  0.5  0.5
[6,] 0.00 0.00 0.00 0.00 0.33  0.0  0.5
[7,] 0.00 0.00 0.00 0.00 0.33  0.5  0.0

> tmp
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.15 0.15 0.15 0.15 0.15 0.15 0.15
[2,] 0.20 0.20 0.20 0.20 0.20 0.20 0.20
[3,] 0.15 0.15 0.15 0.15 0.15 0.15 0.15
[4,] 0.15 0.15 0.15 0.15 0.15 0.15 0.15
[5,] 0.15 0.15 0.15 0.15 0.15 0.15 0.15
[6,] 0.10 0.10 0.10 0.10 0.10 0.10 0.10
[7,] 0.10 0.10 0.10 0.10 0.10 0.10 0.10

tmpのようになってしまうため、どこがグループなんだかさっぱりわからなくなってしまう。

なので一工夫して、最初の仮定、無限回遷移するとグループ内で固まるという状況を、無理矢理つくり上げるのが、MCLである。

The Markov Cluster Algorithm

アルゴリズムは、以下のようになる。


  1. 無向グラフをインプットする。
  2. 累乗のパラメータe、増幅パラメータrを導入する(適当に決める)。
  3. 自己ループさせたければさせる。
  4. 正規化する
  5. e回行列の積を取る
  6. パラメータrにしたがって、グラフの偏りを顕著にする(inflate)
  7. 収束するまで5と6をループする
  8. 結果を解釈してクラスタを決める

6番がキモで、これで無理矢理グラフを分割しやすくしている。

とりあえず実装したので、結果を見る。

#1.行列を初期化
data <- matrix(0,7,7)
data[1,] <- c(0,1,1,1,0,0,0)
data[2,] <- c(1,0,1,1,1,0,0)
data[3,] <- c(1,1,0,1,0,0,0)
data[4,] <- c(1,1,1,0,0,0,0)
data[5,] <- c(0,1,0,0,0,1,1)
data[6,] <- c(0,0,0,0,1,0,1)
data[7,] <- c(0,0,0,0,1,1,0)

#dataと同じ大きさの単位行列を用意する
tmp <- matrix(0,dim(data)[1],dim(data)[2])
diag(tmp) <- 1

#2.パラメータを二つ定める
##行列を何乗するのか
e <- 2
##どのくらいグループを極端に分けるのか
r <- 2

#3.自己ループさせる(やらなくても良い)
diag(data) <- 1

#4.行列を正規化する
data <- t(t(data)*(1/colSums(data)))
FLAG <- data

#7.処理をループさせる
while(TRUE){

  #5.行列をe乗する
  #元論文が見れないので何をe乗するのか定かでない
  #FLAGかdata(FLAG説が濃厚)
  for(i in 1:e){
    tmp<-tmp%*%FLAG
  }

  #6.パラメータrによってグループを顕著にする
  tmp <- tmp^r
  #条件は満たすように正規化
  tmp <- t(t(tmp)*(1/colSums(tmp)))

  #終了判定
  if(sum(FLAG-tmp) < 0.01){
    break
  }
  FLAG <- tmp
}

#8.画像を出力する(デフォ関数じゃないです)
#image(tmp)
image_text(round(tmp,2))

図は間違えてdataをインフレーションさせてしまったものである。

本当はもっと顕著に偏りが生まれるけれど、撮り直すのがめんどいのでそのまま。

e<-2
r<-2
f:id:rishida:20130606205932p:plain:w300

e<-2
r<-10
f:id:rishida:20130606210229p:plain:w300

e<-10
r<-2
f:id:rishida:20130606210322p:plain:w300

e<-10
r<-10
f:id:rishida:20130606210418p:plain:w300

見れば分かるが、1,2,3,4のグループと、5,6,7のグループが出来ている。

アルゴリズムの6番は、全部の値をそれぞれr乗することによって、差を顕著にしているのである。

> v
 [1] 0 1 2 3 4 5 6 7 8 9
> v/sum(v)
 [1] 0.00000000 0.02222222 0.04444444 0.06666667 0.08888889 0.11111111
 [7] 0.13333333 0.15555556 0.17777778 0.20000000
> v^2
 [1]  0  1  4  9 16 25 36 49 64 81
> v^2/sum(v^2)
 [1] 0.000000000 0.003508772 0.014035088 0.031578947 0.056140351 0.087719298
 [7] 0.126315789 0.171929825 0.224561404 0.284210526

こういうイメージ。