各種のランキングをRでためしてみよう(3)—PH Alogrithm

前回、前々回にひきつづき、先生の論文読んで同じ実験をRでやってます:

Hofuku, Oshima "Measures to Represent the Proparties of Nodes in a Drected Graph"(まだ掲載されてないがacceptはされたそうだ).

この論文であつかっているもの:

この論文ではこのこのように何かアルゴリズム(行列)を考えたとき、どのような行列であれば固有ベクトル＝ランキングベクトルを、べき乗法で一意に決定できるかの必要条件も確認してて、その際のキーワードは「non-negative, primitive and irreducible」みたいです。

ちなみに原始行列(primitive matrix)かどうかってのは、

その行列を隣接行列としてみてグラフをつくり、
あるノードiから出てノードiに戻るためのエッジの数を周期と呼び,
グラフ全体についての周期の集合が、最小の周期の倍数である
(たとえば対角成分が1つでも1ならprimitive)

ということらしい。わかんなくてさっき聞いたんだけど。

PHアルゴリズム

PH=PageRank+HITSみたいな

PHのモデルは、PageRank的な文脈でいえば

「よいページと、にている(共通のリンク元や、共通のリンク元のリンク元が多い)ページは、よいページである。」

という感じだろうか。

あと、いままですっかり無視してきたけど、原始的なPageRankでは、スコアの再分配は一律で、ひとつのページ(ノード)からk個のリンク(エッジ)が張ってある場合は、1/kで分配するというものである。
これに対して、PHではエッジに重みがかかるのでそうとは限らない。
ただ、いずれにしても、ベクトルr(t)の要素 $r(t)_i$ を、r(t+1)の各要素に再分配する際に総量が増減しないように、行列の要素を
$r(t)_{i}v_{i1}+r(t)_{i}v_{i2}+\cdots+r(t)_{i}v_{iN}=r(t)_i$
となるように調整しておく必要があって、具体的には
$v_{i1}+\cdots+v_{iN}=1$
となるように設定するということ。これは下記STEP3&STEP4の、行列の各行をその合計(l1-ノルム)が1になるように正規化する操作にあたる。

(ソース全体は一番下に。)

隣接行列

前回と同じ。

> n <- t(cbind(
 + c(0,1,0,0,0,0),
 + c(0,0,1,1,0,0),
 + c(0,0,0,0,0,0),
 + c(0,0,1,0,0,1),
 + c(0,0,0,1,0,0),
 + c(0,0,0,0,0,0)))
> 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

STEP1(add paths to nodes reached within 2-step)

まずPHアルゴリズムでは、nをそのまま使うのではなく、
$M=N+kN^2$
みたいな行列をつくる。
各項のバランスをとるパラメータ(k)がついているのだけど、ここではこれを1として無視する。すると、これは単純に、

それぞれの丸から2ステップで行ける丸に、すべて直接の矢印を追加してしまう

操作に相当する。それがいくつもあった場合はその合計本数になる。

> m<-n+n%*%n
> library(igraph)
> plot(graph.adjacency(m)) # 毎度思うけどなぜRのくせにゼロオリジンなの…

> m
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    1    1    1    0    0
[2,]    0    0    2    1    0    1
[3,]    0    0    0    0    0    0
[4,]    0    0    1    0    0    1
[5,]    0    0    1    1    0    1
[6,]    0    0    0    0    0    0

STEP2(authority matrix)

さてMから直接ランキングを求めるのではなく
$U=^TMM$
として行列Uをつくる(HITSのauthority行列っぽい)。
u[i,j](=u[i,j])は「ノードiとノードjが共通のリンク元(およびリンク元のリンク元)をいくつもつか」を表す。

> u<-t(m)%*%m
> plot(graph.adjacency(u))

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

STEP3&STEP4(normalizing of rows)

つぎに、全要素ゼロの行を1で埋めてから、これの各行を正規化する。
(PageRankっぽい。)

> v<-u
> v[rowSums(v)==0,]<-1
> v<-v/rowSums(v)
> 
> vv<-0.1*(1/6)+0.9*v
> w<-t(vv)
> v
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
[1,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
[2,] 0.0000000 0.3333333 0.3333333 0.3333333 0.0000000 0.0000000
[3,] 0.0000000 0.0625000 0.4375000 0.2500000 0.0000000 0.2500000
[4,] 0.0000000 0.1000000 0.4000000 0.3000000 0.0000000 0.2000000
[5,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
[6,] 0.0000000 0.0000000 0.4444444 0.2222222 0.0000000 0.3333333
> w
          [,1]       [,2]       [,3]       [,4]      [,5]       [,6]
[1,] 0.1666667 0.01666667 0.01666667 0.01666667 0.1666667 0.01666667
[2,] 0.1666667 0.31666667 0.07291667 0.10666667 0.1666667 0.01666667
[3,] 0.1666667 0.31666667 0.41041667 0.37666667 0.1666667 0.41666667
[4,] 0.1666667 0.31666667 0.24166667 0.28666667 0.1666667 0.21666667
[5,] 0.1666667 0.01666667 0.01666667 0.01666667 0.1666667 0.01666667
[6,] 0.1666667 0.01666667 0.24166667 0.19666667 0.1666667 0.31666667

STEP2(power method):

rを初期化。(r <- sign(1:m)/mとかも考えたが、1行で書けるもっとすてきな方法もあるんじゃないかな！)

> r<-vector()
> r[1:6]<-1/6

for文でべき乗法の繰り返しを計算する。rがほとんど変化しなくなる(収束する)までやる。

> r<-vector()
> r[1:6]<-1/6
> j<-1
> re<-vector()
> re[j]<-1
> while(0.000000001 r
           [,1]
[1,] 0.02380952
[2,] 0.09704873
[3,] 0.38262023
[4,] 0.25113558
[5,] 0.02380952
[6,] 0.22157640
> order(r,decreasing=TRUE)
[1] 3 4 6 2 1 5

のように一応、それっぽいスコアとランキングが求まった(1と5はタイ)。20回目からは一応、表示桁では変化なし(下位の桁はまだ動いている)。

このグラフだと、元祖PageRankおよびHITS(authority)のところでやってみた結果と同じですね。:

> order(r,decreasing=TRUE)
[1] 3 4 6 2 1 5

ところで、「r==p%*%r」がすべてTRUEになるまで、完全な収束する判定もできそうですが、そもそも完全に一致するのは理論的に無限大の繰り返し後だし、なかなか(あるいはまったく)止まらなくなります。
こういう風になにかベクトル変化量を定義すると、単調減少しますから、それがしきいちを下回るまで回すとかすればいいんじゃないですかね。(さらに、規定回数を超えても収束しないなら、breakする文も書くべきかもしれませんが。)
上のコードでは、ためしに変化量の合計をベクトルに蓄積してみたのでplotしてみます。plot()するには、21回で止まったので:

> plot(1:21,re[1:21],"l")
> plot(1:24,log(re),"o")

です。(縦がlogのplotもついでに。直線になってるんで、精度を求めると、指数関数的に繰り返し数が増えるということかな。)

参照してる論文には、これのあとに、ノードi-ノードjの、対応するベクトルv[,i]とv[,j]のcos距離を計算して行列にして…とかが書いてあるけどとりあえずここまで。

で、こういうランキングをして、上位からいくつかのノードを選び、それぞれを中心にしたクラスタリングとかを、僕はやるそうです。
一応、べき乗法のところが差分の式にできるから、極限にとばして、連続関数みたいに偏微分しちゃって、制御(工学)みたいな安定がどうのこうの、的なことできないっすかね。とか、いうだけいいましたが、そもそも、おれ、微積できないので微分方程式なんかとんでもないのであった。

Rのソース

library(igraph)

m<-n+n%*%n
plot(graph.adjacency(m))

u<-t(m)%*%m
plot(graph.adjacency(u))

v<-u
v[rowSums(v)==0,]<-1
v<-v/rowSums(v)

vv<-0.1*(1/6)+0.9*v
w<-t(vv)

r<-vector()
r[1:6]<-1/6
j<-1
re<-vector()
re[j]<-1
while(0.000000001