各種のランキングをRでためしてみよう(2)—HITS(authority)

きのうの続きを。
それにしても、ほんとに行列積、行列ベクトル積くらいしか使ってないんですよ、これ。
偉大ですね！　微積とかわかんねーもんな、おれ。

一連のエントリは、先生の論文読んで、同じ実験をRでやっています:

Hofuku, Oshima "Measures to Represent the Proparties of Nodes in a Drected Graph"(未掲載？).

この論文であつかうもの:

原始的なPageRank: きのうの書いた
HITS(using H1 called "authority matrix"): 今回
"PH Algorithm"(上記論文で提案)

HITS(H1を使うやつ=authority)

概要

HITSやります。
HITSというのは、「authority」と「hub」という概念を持ち出して、authorityとしてみたときのランキングと、hubとしてみたときのランキングをそれぞれ計算します。
ほんとはグラフ理論のことばで書くべきかと思って、前エントリは妙な言葉遣いになってましたが、あきらめた。
グラフを図で示す場合、丸はウェブページ、矢印はあるページから別のページへのリンクということにする。

Rで書いたソース全体は末尾に示したとおり。

HITSのモデルは

「よいhubページからリンクされているページはよいauthorityページである。」
および「よいauthorityページにリンクを張っているページはよいhubである。」

という仮説に従っている。

また、PageRankとちがって、ほんとは、クエリに関連するページの集合内で求めるべきなのだけど(クエリに出現した単語をひとつでも含むページをあつめた集合とか)、ここではその作業ははぶいて、ページの集合はあらかじめ与えられているものとしている。

ここで「HITS」「HITS」っていってますが、「PageRank」同様に、やっぱいろいろあるのかもしれない(もう、出てからけっこう経つし。たとえば、たしか、便所先生(@y_benjo)がはてダでco-HITSという言葉をつかってた気がして、気になるのだが)。
ただ、ここでみている論文では、とくに歴史的経緯などについて注釈がついてないので、おそらくこれがいちばん古い(単純な)やつじゃないかと思う。

隣接行列

前回と同じ。
ただし、ちょっと趣向を変えて、c()で列ベクトル作っておいてcbind()する方法をつかってみた。

> 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
> library(igraph)
> plot(graph.adjacency(n))

(どうでもいいんだけど、なんか最後の行、graphオブジェクトをその場で生成しているせいか、plot()の結果が毎度かわりますね…。)

STEP1(行列h1とh2):

authorityスコア算出につかう行列(authority行列という)h1と、hubスコア算出につかう行列(hub行列という)h2を、それぞれ作ります。
違いは「転置n・n」か、「n・転置n」か。(転置というのは…対角成分(左上-右下)のラインに対して対称に、右上と左下をくるっと裏返すような感じの操作ですね、はい)。

> h1<-t(n)%*%n
> h1
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    0    0    0    0    0
[2,]    0    1    0    0    0    0
[3,]    0    0    2    1    0    1
[4,]    0    0    1    2    0    0
[5,]    0    0    0    0    0    0
[6,]    0    0    1    0    0    1
> h2<-n%*%t(n)
> h2
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    0    0    0    0
[2,]    0    2    0    1    1    0
[3,]    0    0    0    0    0    0
[4,]    0    1    0    2    0    0
[5,]    0    1    0    0    1    0
[6,]    0    0    0    0    0    0

ただし、Rにおいて「%*%」は行列-行列積。
t(n)*nとすると、Rは同じ場所の要素同士の積を計算してくれちゃいますので、注意が必要(まちがえてゼロ行列になったり…した)。

STEP2(h1を使ってべき乗法):

rを初期化。

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

h1の規格化をしてないので、少なくとも、「r<-h1%*%r」を繰り返すと単調増加してしまっていくつかの要素は無限大に発散することは予想できる。そもそも、どうも以下のfor文の繰り返しでも、初期値に依存せず一点に収束するかの保証はないらしいのだが。
そこで、
$\bf{r}_{(t+1)}:=\frac{\bf{r}_{(t)}}{\sum^n{r_{(t)i}}$
として、毎回、前要素の合計が1になるようにしてみた。
こうすると、

> r<-vector()
> r[1:6]<-1/6
> for(i in 1:6)r<-h1%*%(r/sum(r))
> r/sum(r)
            [,1]
[1,] 0.000000000
[2,] 0.000311042
[3,] 0.445101089
[4,] 0.356143079
[5,] 0.000000000
[6,] 0.198444790
> order(r,decreasing=TRUE)
[1] 3 4 6 2 1 5

のように一応、それっぽいauthorityとしてのスコアとランキングが求まった(1と5はタイ)。
同様にh2を使って繰り返し計算をするとhubとしてのスコアスコアとランキングも求められる。
rは絶対値最大の固有ベクトルと、少なくとも平行なものになっている…といいなあ…。

そんなときはこれを:
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にも言及してる。とても網羅的で便利。教科書みたいだ。

いま流し読みしたら

the diﬀerent limiting authority (and hub) vectors can be produced by diﬀerent choices of the initial vector
(p.140)

つまり、ランキングは、繰り返しにつかうベクトルの初期値に依存することがある、とか書いてあったけど気にしない！
あときのうひどい説明を書いたペロン-フロベニウスの定理についても、引用しておきます:

The Perron–Frobenius theorem [46] ensures that an irreducible nonnegative matrix possesses a unique normalized positive dominant eigenvector, called the Perron vector.
(p.141)

HITSだとこれの条件を満たさないことがあるということです。

とりあえず、この計算では、authorityランキングと前回のPageRankのランキング(下記)と一致した。

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

余談: 行列h1、h2について

h1、h2の各要素がなにを意味しているのか。

h1=t(n)*nはnの列ベクトル同士の積である。要素n1[i,j]はページiとページjとが、いくつ同一のリンクを含むか(リンク先ページk(k=1...6))を表す。(nを縦にみて同じ場所に1が立っている数をかぞえるのとおなじ。)
h2=n*t(n)はnの行ベクトル同士の積である。要素n1[i,j]はページiとページjとが、いくつ共通のリンク元ページk(k=1...6)をもっているかを表す。(nを横にみて同じ場所に1が立っている数をかぞえるのとおなじ。)

ちなみにh1、h2を隣接行列とみて

> plot(graph.adjacency(h1,"undirected"))

> plot(graph.adjacency(h2,"undirected"))

みたいなplot()もできる。
ただし、h1、h2はどちらも同じサイズの対称行列なので、無向グラフにした。
この図の線(エッジ、辺)は、authorityとして(被リンク基準)の近さ、hubとして(リンク先基準)の近さをそれぞれ表していて、もとのウェブページの直接のリンク-被リンク関係(隣接行列n)の図とは関係ない。
(隣接行列h1、h2には、それぞれエッジに重みがついている(0、1以外の要素がある)から、うまくやれば線の太さとかに置き換えて図にできるはず。ちょっとplot()のこつがわからないのでできてないんだけど。これオフラインで書いてるので、ググれないんです。)

Rのソース

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)))

library(igraph)
plot(graph.adjacency(n))

h1<-t(n)%*%n
h2<-n%*%t(n)

r<-vector()
r[1:6]<-1/6
for(i in 1:6)r<-h1%*%(r/sum(r))
r/sum(r)

(あーつかれた…。もう1つ残ってんだけど、どうしようかなー…。)