前エントリにひきつづき、研究のための体制立て直しの様子を、垂れ流していきたいと思います!
隣接行列
こんどは隣接行列 (adjacency matrix) からグラフをつくる。前エントリでみたRion778さんのページを、ひきつづき参考にします。
graph.adjacency()関数により,隣接行列 (adjacency matrices) からグラフオブジェクトを作成できる.対角成分が0でなければループになる.
そこでまず、
> n <- c( + 0,1,0,0,0,0, + 0,0,1,1,0,0, + 0,0,0,0,0,0, + 0,0,1,0,0,1, + 0,0,0,1,0,0, + 0,0,0,0,0,0) > n <- t(matrix(n, nrow=6)) > n [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0 1 0 0 0 0 [2,] 0 0 1 1 0 0 [3,] 0 0 0 0 0 0 [4,] 0 0 1 0 0 1 [5,] 0 0 0 1 0 0 [6,] 0 0 0 0 0 0
という行列nを手入力した(ほんとはCSVファイルなどからも読めたはずなので、実際にはこのようなやりかたは、何度もデータを換えて実験するときには、あまりおすすめしないが)。
graph.adjacency {igraph} R Documentation
Create graphs from adjacency matrices...
Usage
graph.adjacency(adjmatrix, mode=c("directed", "undirected", "max",
"min", "upper", "lower", "plus"), weighted=NULL, diag=TRUE,
add.colnames=NULL, add.rownames=NA)
> plot(graph.adjacency(n))
ここで「あ、ちげえ。plotして気づいたけど、nが隣接行列の転置になってるじゃん!」……などと、いうこともある(直した)。
Hofuku, Oshima "Measures to Represent the Proparties of Nodes in a Drected Graph"(あれ未発表?そういえば、これ書いてだいじょぶか?).
この論文のあつかうもの:
以降、この論文と同じ実験をRでやって同じ結果が得られるかをみます。
とりあえず、きょうはPageRankのかんたんなやつを。
追記:
いまさら「RでPageRank」とかググったら、igraph内に、そのものすばりgraph.rankというものが…。
- Rでページランクの計算 - arupaka-_-arupakaの日記
- http://rbloger.blog51.fc2.com/blog-entry-17.html" title=統計ソフトRのブログ PageRank
Calculates the Google PageRank for the specified vertices.
Usage
page.rank (graph, vids = V(graph), directed = TRUE, damping = 0.85,
weights = NULL, options = igraph.arpack.default)
page.rank.old (graph, vids = V(graph), directed = TRUE, niter = 1000,
eps = 0.001, damping = 0.85, old = FALSE)
ぎゃふん!
原始的なPageRank
元祖のPageRankやります。
このPageRankは「よいページからリンクされているページはよいページである。」みたいに、リンクを単純な相互評価としてとらえる仮説に基づいて、みんなの意見を統一したような評価を得るモデルです。訪問者数のある種の期待値ともいえるそうな。
参照してる論文の書き方、どうも少し特殊な気がするので、少しはしょって書くけど。
ソース全体はこんな感じ。
n <- c( 0,1,0,0,0,0, 0,0,1,1,0,0, 0,0,0,0,0,0, 0,0,1,0,0,1, 0,0,0,1,0,0, 0,0,0,0,0,0) n <- t(matrix(n, nrow=6)) p<-n p[rowSums(p)==0,]<-1 p<-p/rowSums(p) pp<-0.1*(1/6)+0.9*p q<-t(pp) r<-vector("numeric",6) r[]<-1/6 for(i in 1:10)r<<-q%*%r r
以下、詳細。
STEP1&STEP2(規格化):
各行のl1-ノルム(この場合、単純な合計)に関して、1になるように調整する(各行を、おのおののl1-ノルムで割る)(STEP1)。
ただし、いきなりn/rowSums(n)とやりたくなるが、これだとゼロ割りしてやっかいなので、合計がゼロのとこだけ1で割るようにしておきます(3行目)ぜんぶゼロの行はまず非ゼロの、適当な値を一律に入れてしまう(STEP2)。
ちなみに、以前、参考にした論文も:
- Amy N. Langville, Carl D. Meyer 'A Survey of Eigenvector Methods for Web Information Retrieval', SIAM REVIEW c Vol. 47, No. 1, pp. 135–161 (Society for Industrial and Applied Mathematics, 2005).
- http://meyer.math.ncsu.edu/Meyer/PS_Files/Survey.pdf
- HITS、PageRank、SALSAとか「絶対値最大の固有値に対する固有ベクトル」をランキング指標として使う方面の手法の、まとめ。べき乗法、ペロン=フロベニウスの定理とか、そんな感じ。
- 同じ著者による本の日本語版もある
Google PageRankの数理 ―最強検索エンジンのランキング手法を求めて―
- 作者: Amy N.Langville,Carl D.Meyer,岩野和生,黒川利明,黒川洋
- 出版社/メーカー: 共立出版
- 発売日: 2009/10/10
- メディア: 単行本
- 購入: 4人 クリック: 249回
- この商品を含むブログ (23件) を見る
これによると、ゼロの行にいれるのは、「 、ただしnはPのオーダー」だそうな。具体的にはよくわかんないので1/n(=1/6)にしておく。これをpとする(よく考えたら、各行の合計が1にならないとまずいのです)。
> p<-n > tmp<-rowSums(n) > p[tmp==0,]<-1 > tmp[tmp==0]<-6 > p<-p/tmp > p [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [2,] 0.0000000 0.0000000 0.5000000 0.5000000 0.0000000 0.0000000 [3,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 [4,] 0.0000000 0.0000000 0.5000000 0.0000000 0.0000000 0.5000000 [5,] 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 [6,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
STEP3(べき乗法を収束させる条件):
行列pをともかく既約にする(グラフを強連結にする)ために、
みたいなことをする。ただし、Eは全要素が1で、c=0.9とする。Rでは右辺をp、左辺はppと書いている。
(いちおう、より一般的にであるEを用いて、vを動かしてもっと細かく制御する方法もあるらしい。けど、めんどいのでやりません。)
> pp<-0.1*(1/6)+0.9*p > pp [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.01666667 0.91666667 0.01666667 0.01666667 0.01666667 0.01666667 [2,] 0.01666667 0.01666667 0.46666667 0.46666667 0.01666667 0.01666667 [3,] 0.04166667 0.04166667 0.04166667 0.04166667 0.04166667 0.04166667 [4,] 0.01666667 0.01666667 0.46666667 0.01666667 0.01666667 0.46666667 [5,] 0.01666667 0.01666667 0.01666667 0.91666667 0.01666667 0.01666667 [6,] 0.04166667 0.04166667 0.04166667 0.04166667 0.04166667 0.04166667 > q<-t(pp)
STEP4(べき乗法で固有ベクトル):
> r<-vector("numeric",6) > r[]<-1/6 > r<-q%*%r > r [,1] [1,] 0.06666667 [2,] 0.21666667 [3,] 0.21666667 [4,] 0.29166667 [5,] 0.06666667 [6,] 0.14166667 > r<-q%*%r > r [,1] [1,] 0.07041667 [2,] 0.13041667 [3,] 0.29916667 [4,] 0.22791667 [5,] 0.07041667 [6,] 0.20166667 > r<-q%*%r > r [,1] [1,] 0.09179167 [2,] 0.15516667 [3,] 0.25304167 [4,] 0.21385417 [5,] 0.09179167 [6,] 0.19435417 > r<-q%*%r > r [,1] [1,] 0.08377604 [2,] 0.16638854 [3,] 0.24983542 [4,] 0.23621354 [5,] 0.08377604 [6,] 0.18001042 > r<-q%*%r > r [,1] [1,] 0.08114354 [2,] 0.15654198 [3,] 0.26231448 [4,] 0.23141682 [5,] 0.08114354 [6,] 0.18743964 > r<-q%*%r > r [,1] [1,] 0.08412978 [2,] 0.15715897 [3,] 0.25871124 [4,] 0.22760286 [5,] 0.08412978 [6,] 0.18826735 > r<-q%*%r > r [,1] [1,] 0.08371346 [2,] 0.15943026 [3,] 0.25685628 [4,] 0.23015180 [5,] 0.08371346 [6,] 0.18613474
と、まあこんな感じか。7回ですか。
左を一位として、ノードの番号を表示してみると
> order(r,decreasing=TRUE) [1] 3 4 6 2 1 5
となる(1と5はタイ)。
さっきも見た画像ですが、これに「重要なノードからのエッジがあるノードは重要である」というモデルを適用すると、こういう順番になるわけです。
ちなみに、何度「r<<-q%*%r」を繰り返しても
> sum(r) [1] 1
になってるはずで、実際
「0.08371346+0.15943026+0.25685628+0.23015180+0.08371346+0.18613474=1」です。
これははじめに行列の各行のl1-ノルム(合計)が1になるように規格化したからです。また「絶対値最大の固有値が1だ」ということとも関連している。上記の無限の繰り返しによって得られたものが、そのqの固有値1に対する固有ベクトル。
あと、増える要素は単調に増えるし、減る方も同様(たぶん)。なので、エネルギー保存則みたいなものというか、そういう感じですたぶん。要素がすべて正の実数の既約な行列の、絶対値最大の固有値がべき乗法で求まることのの保証とかは、ペロンさんとフロベニウスさんにお聞きください。
余談
で、いまさらだけど、じつはRにはeigenという関数があって、行列のサイズが小さければ
> eigen(cbind(c(1,-1),c(-1,1))) $values [1] 2 0 $vectors [,1] [,2] [1,] -0.7071068 -0.7071068 [2,] 0.7071068 -0.7071068
みたいに一行で求まります。べき乗法の必要がない、ふつうの場合はたぶんこちらを使います。
とかいって、でも固有ベクトルが一致しないんですよね、なぜか…。
> eigen(q) $values [1] 1.000000e+00+0.000000e+00i -1.789082e-01+4.677976e-01i [3] -1.789082e-01-4.677976e-01i -2.421836e-01+0.000000e+00i [5] 3.801593e-09+0.000000e+00i -3.801593e-09+0.000000e+00i $vectors [,1] [,2] [,3] [1,] -0.1896813+0i -0.1616564-0.1997846i -0.1616564+0.1997846i [2,] -0.3603945+0i -0.3932094+0.1997846i -0.3932094-0.1997846i [3,] -0.5870163+0i 0.5548659+0.0000000i 0.5548659+0.0000000i [4,] -0.5225721+0i -0.0993472+0.4656472i -0.0993472-0.4656472i [5,] -0.1896813+0i -0.1616564-0.1997846i -0.1616564+0.1997846i [6,] -0.4248388+0i 0.2610036-0.2658626i 0.2610036+0.2658626i [,4] [,5] [,6] [1,] -0.1861799+0i 2.515346e-16+0i -2.355139e-16+0i [2,] 0.5056998+0i 1.194726e-08+0i 1.194725e-08+0i [3,] -0.3195199+0i 7.071068e-01+0i -7.071068e-01+0i [4,] -0.4339382+0i -5.973627e-09+0i -5.973629e-09+0i [5,] -0.1861799+0i -5.973627e-09+0i -5.973628e-09+0i [6,] 0.6201181+0i -7.071068e-01+0i 7.071068e-01+0i