グラフ(ネットワーク状)構造のクラスタリング(グループ分け)基礎:The Markov Cluster Algorithm(MCL)まとめた
an introduction to the MCL algorithm
めちゃくちゃ基本らしい。けど微妙に関係ないので、さわりだけ読む。
無意味にグラフ理論の基礎に触れつつ、エッジリストから高速にクラスタを組む話をする。
元論文は有料だったので諦めた。ので、何か微妙に間違っているかもしれない。コードあるけどR言語です。
グラフ構造とは
グラフ構造とは、網目状の構造である。例えば、鉄道とか道路とかネットとか神経回路がそれに相当する。
例えば、道路の例を考えよう。
丸(ノード)が交差点で、線(エッジ)が道路になる。
この図は、「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と呼ぶ。
は、jから出発して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で割れば良い。
したがって、列和で列を割ることにする。
この結果として、には、jからiに移動する確率が入る。
この確率による移動の連鎖が、一次マルコフ連鎖という代物である。
最初に各交差点にいる確率(今回は単位行列を暗に仮定)と、iにいる人が、jに移動する確率は定まっている(マルコフ連鎖)
そしてその移動確率は、今いる場所にだけ依存する(一次)
グラフにおけるグループ
さて、は、jからiに移動する確率であるとした。
ここで、行列X同士の積を取る。
とかすると、
つまり、は、jからkに行ってkからiに行く確率になるので、結局は、jから出発してなんだかんだあってiに到着する確率になる。
今回の場合は消滅したりしないので、jから出発した人はどこかしらにいる。
つまり、のあるjに関して、全てのiについて足し合わせると(行列の列和を取ると)、総和は1になる。
もも列和は常に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がもうひとつのグループに見える。
このグループを単純な数値計算で簡単に見つけるのが、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
アルゴリズムは、以下のようになる。
- 無向グラフをインプットする。
- 累乗のパラメータe、増幅パラメータrを導入する(適当に決める)。
- 自己ループさせたければさせる。
- 正規化する
- e回行列の積を取る
- パラメータrにしたがって、グラフの偏りを顕著にする(inflate)
- 収束するまで5と6をループする
- 結果を解釈してクラスタを決める
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
e<-2
r<-10
e<-10
r<-2
e<-10
r<-10
見れば分かるが、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
こういうイメージ。