グラフ・ネットワーク分析で遊ぶ(2)：最短経路長など

前回の記事に引き続き主に{igraph}の各関数で遊びながらグラフ理論・ネットワーク分析を学ぶこのシリーズですが、今回は様々なノード間の特徴量について見てみます。もちろん今回も参考文献はこちら。

ネットワーク分析 (Rで学ぶデータサイエンス 8)

データセットは前回適当に生成したグラフのものと、C elegansと、さらに以前使った『レ・ミゼラブル』の人物相関図を対比のために併用しようと思います。

最短経路長とダイクストラ法

前回適当に生成したグラフを今回も使ってみましょう。

> d
      [,1] [,2]
 [1,]    1    8
 [2,]    3    1
 [3,]    1    2
 [4,]    6    2
 [5,]    3    8
 [6,]    5    3
 [7,]    6    3
 [8,]    8    3
 [9,]    3    4
[10,]    8    4
[11,]    2    5
[12,]    4    5
[13,]    8    5
[14,]    1    6
[15,]    2    6
[16,]    4    6
[17,]    5    6
[18,]    4    7
[19,]    6    7
> library(igraph)
> g0<-graph.edgelist(as.matrix(d),directed=T)
> set.seed(71)
> plot(g0,layout=layout.fruchterman.reingold)
# パッケージをアップデートしたので配色が前回と違う

こんな感じになります。で、最初のお題は「最短経路長」(shortest path length)。上のプロットを見れば一目瞭然だと思いますが、例えばノード1から2への最短経路長はどう見ても1ですね（1エッジで到達する）。一方で、ノード2から8は3もあります（3エッジかかる）。有向グラフだと当然ながら迂回するような形でしかあるノードから別のノードへと到達できないことがあるので、グラフが大規模になった場合のその計算はちょっと厄介です。

そこで、最短経路長の算出法としてよく用いられるのがダイクストラ法。Wikipedia日本語版記事ではこんなアナロジーで説明されています。

簡単の為、G=(V,E)の各頂点v,w∈Vに対し、vとwを結ぶ辺は最大一つしかないとする。 vとwを結ぶ辺があれば、その辺を(v,w)と書く。

最短経路問題は、ビー玉と紐を用いて解くことが出来る。 まずビー玉を頂点、紐を辺にするグラフを工作する。 グラフを板の上に置き、スタートの頂点にあたるビー玉だけをつまむ。 グラフが置かれている板を取り除くと、グラフは自由落下を始めるが、 スタートにあたるビー玉を持っているので、スタート地点から近いビー玉から順に落下が止まる。 ゴールにあたるビー玉が止まったとき、ゴールにあたるビー玉はスタートにあたるビー玉まで紐で一直線で結ばれている。 この直線が最短経路である。

落下が止まった頂点vに対し、vを支えている頂点wが存在する。 wをvの「直前の頂点」と呼ぶことにする。 以下簡単のため、直前の頂点は1つしかないと仮定して話を進めるが、 一般の場合も同様である。

ダイクストラのアルゴリズムは、上述の解法をコンピュータ上でシミュレートしたものである。 グラフG=(V,E)、スタートとなる頂点s、および各辺eの長さlength(e)が入力として与えられると、 このアルゴリズムは上述の解法をシミュレートし、「落下」が停止した頂点から順に、その頂点の直前の頂点が何であるかを出力する。 ゴールとなる頂点の「落下」が停止したところで、ダイクストラのアルゴリズムを止めれば、 あとはゴールから順に、直前の頂点、さらにその直前の頂点とたどっていくことで、 ゴールとスタートをつなぐ最短経路(の一つ)を求めることができる。

『ネットワーク分析』pp.22-25にはそのアルゴリズムを追ったRスクリプトが載っているので、そのまま写経して実行してみます。

> dijkstra<-function(matrix,vertex){
+     L<-matrix
+     L[which(L==0)]<-Inf
+     diag(L)<-0
+     n<-nrow(matrix)
+     d<-rep(Inf,n)
+     d[vertex]<-0
+     M<-1:n
+     M<-M[-vertex]
+     i<-vertex
+     while(length(M)>0){
+         for (j in 1:n)
+             d[j]<-min(d[j],d[i]+L[i,j])
+         i<-M[which(d[M]==min(d[M]))[1]]
+         M<-M[-which(M==i)]
+     }
+     d
+ }
> A0<-get.adjacency(g0)
> dijkstra(as.matrix(A0),1) # ノード1からの距離
[1] 0 1 2 2 2 1 2 1
> dijkstra(as.matrix(A0),8) # ノード8からの距離
[1] 2 3 1 1 1 2 2 0

プロットで見た通りの結果になっていることが分かります。ダイクストラ法はエッジの重み付けが非負の時にのみ使える（この例では全部重みを等しく全て1としてますが）方法で、その他の場合にはそれぞれ対応するアルゴリズムがあります。{igraph}のshortest.paths関数はそれぞれの条件に合わせて最適のアルゴリズムを自動的に選択し、最短経路長の算出に用います。

> gs0<-shortest.paths(g0,mode="out")
# 有向グラフの場合はmode引数はoutを指定する
> gs0
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    0    1    2    2    2    1    2    1
[2,]    3    0    2    3    1    1    2    3
[3,]    1    2    0    1    2    2    2    1
[4,]    3    2    2    0    1    1    1    3
[5,]    2    2    1    2    0    1    2    2
[6,]    2    1    1    2    2    0    1    2
[7,]  Inf  Inf  Inf  Inf  Inf  Inf    0  Inf
[8,]    2    3    1    1    1    2    2    0

自分から出て行くエッジを持たないノード7からの最短経路長が全てInfになっている点に注意。また1行目と8行目が上記のdijkstra関数と同じ結果を与えていることも分かります。

最短経路長と可視化

ところで、前回の記事でいくつか力指向グラフ描画アルゴリズムを取り上げましたが、その中で例えばFruchterman - Reingoldアルゴリズムなどは「ノード間の距離に関連するように配置する」という話をしたかと思います。この「距離」はもちろん最短経路長であって、例えばノード1から近くなるほどノードを表す円が大きくなるように描画するとこんな感じになります。

> set.seed(71)
> plot(g0,vertex.size=50/(gs0[1,]+0.5))

多少あいまいですが、概ねプロット上でのノード1までの距離が近いノードほど大きく描画されているのが分かるかと思います。同じことをC elegans、『レ・ミゼラブル』のデータでもやってみましょう。

> g1<-read.graph("celegansneural.gml",format="GML")
> g2<-read.graph("lesmis.gml",format="GML")
> gs1<-shortest.paths(g1,mode="out")
> gs2<-shortest.paths(g2,mode="out")
> set.seed(71)
> plot(g1,vertex.size=50/(gs1[1,]+0.5))
> set.seed(71)
> plot(g2,vertex.size=50/(gs2[1,]+0.5))

ノード1（最も大きな円）からプロット上で見て離れていくほど、概ねどのノードも小さくなっていくのが分かるかと思います。ということで、これでついでにFruchterman - Reingoldアルゴリズムの特性も何となく分かりました、ということで（笑）。