PageRankの計算 DouKaku?

ページランク - Wikipediaを求めるプログラムを作ってください。(PageRankはGoogleの商標です)

詳しい導出方法はGoogle の秘密 - PageRank 徹底解説の3章に載っていますが、簡単に説明すると

ページがn枚ある場合、n×nの2次元配列を用意する。
ページiからページjにリンクがある場合、mat[j][i] = 1 / num_link[i]とする。ただしnum_link[i]はi番目のページから出ているリンクの総数。
行列計算ライブラリを使ってできあがった2次元配列の固有値、固有ベクトルを求める。
出力された固有ベクトルから合計が1になるようなPageRankを算出する

という流れになります。

Pythonで表現すると下のようになります。 ?????の部分は空行を入れて10行でしたので何百行ものコードになってしまった場合はきっとお題の趣旨から外れていると思います。このお題の趣旨は「行列計算ライブラリを使って」PageRankを計算することなので、自力で固有値の計算を実装することは求められていません。

data = {
    1: [2, 3, 4, 5, 7],
    2: [1],
    3: [1, 2],
    4: [2, 3, 5],
    5: [1, 3, 4, 6],
    6: [1, 5],
    7: [5],
}

?????

print pagerank
# [0.303514376997, 0.166134185304, 0.140575079872,
#  0.105431309904, 0.178913738019, 0.0447284345048,
#  0.0607028753994]

このお題はところてんさんの「行列演算系のお題が欲しい」という要望を元に考えたものです。ありがとうございました。

一番乗り？

出力は

{
0.3035143769968051,
0.16613418530351443,
0.1405750798722045,
0.10543130990415332,
0.17891373801916935,
0.04472843450479228,
0.06070287539936102
}

という感じです。

JAMAを使いました。

see: JAMA

import Jama._

def pagerank(data:Map[int, List[int]]):Array[double] = {
  val m = (1 to data.size).map(x=>Array.make(data.size, 0d)).toArray
  for(val i <- data) {
    for(val j <- i._2) {
      m(j-1)(i._1-1) = 1.toDouble / i._2.size
    }
  }

  val m2 = new Matrix((new Matrix(m)).eig.getV.getArray.map(x=>x(0)), data.size)
  val norm1 = m2.norm1
  m2.getArray.map(x => x(0) / norm1)
}
println(pagerank(Map(
  1 -> List(2,3,4,5,7),
  2 -> List(1),
  3 -> List(1,2),
  4 -> List(2,3,5),
  5 -> List(1,3,4,6),
  6 -> List(1,5),
  7 -> List(5)
)).mkString("{\n",",\n","\n}"))

c-wrapper経由でGSLを呼んでみました。
CのAPIそのままなので冗長になってしまう感じですね。
実用的にはこの上に高レベルAPIをかぶせたGauche-gslみたいなのがあると良さそうです。

実行例 (横にはみ出ないように出力は整形してあります):
gosh> (pagerank '((1 . (2 3 4 5 7))
                  (2 . (1))
                  (3 . (1 2))
                  (4 . (2 3 5))
                  (5 . (1 3 4 6))
                  (6 . (1 5))
                  (7 . (5))))
#f64(0.30351437699680506 0.16613418530351437 0.14057507987220447
     0.10543130990415332 0.17891373801916935 0.04472843450479234
     0.06070287539936102)

要Gauche-0.8.11, GSL-1.9です。

(use srfi-1)
(use srfi-42)
(use gauche.uvector)
(use gauche.sequence)
(use c-wrapper)

(c-load '("gsl/gsl_eigen.h" "gsl/gsl_complex_math.h") :libs "-lgslcblas -lgsl")

(define (pagerank data)
  (let* ([dim (apply max (map car data))]
         [dvec (make-f64vector (* dim dim))])
    (do-ec (: r data) (: j (cdr r))
           (set! (ref dvec (+ (* (- j 1) dim) (- (car r) 1))) (/. (length (cdr r)))))
    (let ([eval (gsl_vector_complex_alloc dim)]
          [evec (gsl_matrix_complex_alloc dim dim)]
          [ws   (gsl_eigen_nonsymmv_alloc dim)]
          [mat  (gsl_matrix_view_array dvec dim dim)])
      (gsl_eigen_nonsymmv (ptr (ref mat 'matrix)) eval evec ws)
      (gsl_eigen_nonsymmv_free ws)
      (let* ([maxevi (find-max (iota dim)
                               :key (lambda (k)
                                      (gsl_complex_abs
                                       (gsl_vector_complex_get eval k))))]
             [ev (map-to <f64vector>
                         (lambda (k)
                           (gsl_complex_abs
                            (gsl_vector_complex_get
                             (ptr (ref (gsl_matrix_complex_column evec maxevi) 'vector))
                             k)))
                         (iota dim))])
        (gsl_vector_complex_free eval)
        (gsl_matrix_complex_free evec)
        (f64vector-div ev (f64vector-dot ev (make-f64vector dim 1.0)))))))

標準では、固有値も行列演算もサポートされていないので、自分で計算しています。結果は、合っていると思います。

0.30351438
0.16613419
0.14057508
0.10543131
0.17891374
0.04472843
0.06070287

using System;
class Program
{
  static void Main()
  {
    string data = "2 3 4 5 7,1,1 2,2 3 5,1 3 4 6,1 5,5";
    double[] pageRank = Calc(data);
    foreach (double d in pageRank) Console.WriteLine("{0:F8}", d);
  }
  static double[] Calc(string data)
  {
    string[] ss = data.Split(',');
    double[,] m = new double[ss.Length, ss.Length];
    for (int i = 0; i < ss.Length; ++i)
    {
      string[] tt = ss[i].Split();
      for (int j = 0; j < tt.Length; ++j)
        m[int.Parse(tt[j]) - 1, i] = 1.0 / tt.Length;
    }
    double[] pageRank = new double[ss.Length];
    for (int i = 0; i < pageRank.Length; ++i) pageRank[i] = 1.0 / pageRank.Length;
    double df = 1;
    while (df > 1e-8)
    {
      double[] nxt = new double[ss.Length];
      for (int i = 0; i < ss.Length; ++i)
        for (int j = 0; j < ss.Length; ++j)
          nxt[i] += m[i, j] * pageRank[j];
      df = 0;
      for (int i = 0; i < ss.Length; ++i)
        df = Math.Max(df, Math.Abs(nxt[i] - pageRank[i]));
      pageRank = nxt;
    }
    return pageRank;
  }
}

rb-gslやnarrayのインストールを試みましたが，
エラーが出て面倒になってきた
(ちなみに使っている環境はCygwin)ので
PageRankを求めるアルゴリズムを自分で
実装しました．

一応行列演算のために標準添付されている
Matrixクラスを使っていますが，
論旨に反してごめんなさい．

require "matrix"
require "pp"

def map_to_mat( map )
  size = map.keys.length
  
  con_mat = Array.new( size ) {
    Array.new( size ) { 0.0 }
  }
  
  map.each {| key, val |
    length = val.length
    
    val.each {| i |
      con_mat[ key - 1 ][ i - 1 ] = 1.0 / length
    }
  }
  
  Matrix[*con_mat]
end

def eigen_values( mat, time = 100 )
  time.times {| i |
    mat *= mat
    
    mat_ary = mat.to_a
    mat_ary.map! {| row |
      sum = row.inject( 0.0 ) {| e, t | e + t }
      row.map! {| val | val / sum }
    }
    mat = Matrix[*mat_ary]
  }
  
  mat.to_a[0]
end

## 以下テスト

data = {
  1 => [2, 3, 4, 5, 7],
  2 => [1],
  3 => [1, 2],
  4 => [2, 3, 5],
  5 => [1, 3, 4, 6],
  6 => [1, 5],
  7 => [5],
}

matrix = map_to_mat( data )
pp matrix.to_a

pp eigen_values( matrix )

間違えた…

×論旨
○趣旨

存在する確率の初期ベクトルを等確率と置いて
収束させる方法の方が，#2339より効率的で正攻法ですね．

なんで行列をそのまま回そうと思ったんだろう…

ということで直してみました↓

require "matrix"
require "pp"

def map_to_mat( map )
  size = map.keys.length
  
  con_mat = Array.new( size ) {
    Array.new( size ) { 0.0 }
  }
  
  map.each {| key, val |
    length = val.length
    
    val.each {| i |
      con_mat[ key - 1 ][ i - 1 ] = 1.0 / length
    }
  }
  
  Matrix[*con_mat]
end

def eigen_values( mat, time = 100 )
  size = mat.row_size
  vec = Vector[ *Array.new( size ) { 1.0 / size } ].covector
  
  time.times {| i |
    vec = vec * mat
  }
  
  vec.to_a[0]
end

## 以下テスト

data = {
  1 => [2, 3, 4, 5, 7],
  2 => [1],
  3 => [1, 2],
  4 => [2, 3, 5],
  5 => [1, 3, 4, 6],
  6 => [1, 5],
  7 => [5],
}

matrix = map_to_mat( data )
pp matrix.to_a

pp eigen_values( matrix )

Squeak Smalltalk で。あいにく組み込みの Matrix が貧弱で、固有値や固有ベクトルを求めるメソッドが見あたらなかったので、tell さんの #2339 同様に累乗法で絶対値最大の固有値に対応する固有ベクトルを求めました。あしからず。

| data matrix pageRanks |
data := #((2 3 4 5 7) (1) (1 2) (2 3 5) (1 3 4 6) (1 5) (5)).
matrix := Matrix new: data size element: 0.

data keysAndValuesDo: [:page :refs |
    refs do: [:ref | matrix at: page at: ref put: 1 / refs size]].

pageRanks := [:mat |
    | prev delta |
    prev := mat.
    [   mat := (mat := mat +* mat) / mat sum asFloat.
        ((delta := prev - (prev := mat)) * delta) sum > 1e-12] whileTrue.
    (mat atRow: 1) / (mat atRow: 1) sum].

^pageRanks value: matrix

Javaも標準では行列計算がサポートされていないので、冪乗法で。
固有ベクトルの計算はページランク計算を前提としない、汎用のもののつもり。（ライブラリの代わりなので）
収束しない場合にArithmeticExceptionを投げてるけど、これは手抜き。ホントは（呼び出し元は収束するか事前に判断できないので）チェック例外を投げるべき。

public class PageRank {
    private static final int MAX_REPEAT = 1000;

    public static void main(String[] argv) {
        int[][] data = {{ 2, 3, 4, 5, 7 }, 
                        { 1 }, 
                        { 1, 2 }, 
                        { 2, 3, 5 },
                        { 1, 3, 4, 6 }, 
                        { 1, 5 }, 
                        { 5 }};

        double[][] a = new double[data.length][];
        for (int i = 0; i < a.length; i++) {
            for (int j = 0; j < data[i].length; j++) {
            a[data[i][j] - 1][i] = 1.0 / data[i].length;
            }
        }
        double[] eigVec = calcEigenVector(a);

        double sum = 0.0;
        for (int i = 0; i < eigVec.length; i++) sum += eigVec[i];
        for (int i = 0; i < eigVec.length; i++) System.out.println(eigVec[i] / sum);
    }

    /** べき乗法により、行列の固有ベクトルを求める。 */
    private static double[] calcEigenVector(double[][] mat) {
        final int n = mat.length;
        if(n != mat[0].length) throw new IllegalArgumentException("正方行列じゃないとダメだよ。");
        double eigVec[] = new double[n];
        eigVec[0] = 1.0;

        double eigVal = 1.0;
        
        for(int count=0; count < MAX_REPEAT; count++) {
            double[] newVec = new double[n];
            double df = 0;
            for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) newVec[i] += mat[i][j] * eigVec[j];
            df = Math.max(df, Math.abs(newVec[i] - eigVal * eigVec[i]));
            }
            if(df < 1.0e-8) return eigVec;
            
            double norm = 0;
            for (int i = 0; i < n; ++i) norm += newVec[i] * newVec[i];
            eigVal = Math.sqrt(norm);
            for (int i = 0; i < newVec.length; i++) eigVec[i] = newVec[i] / eigVal;
        }
        throw new ArithmeticException(MAX_REPEAT+"回やっても収束しませんでした。");
    }
}

正直、行列というか数学自体の素養がないので、
お題からリンクされているサイトのコードを、
SciPyを使って内容を理解せずに移植しました。

特に7行めあたりが怪しい。

from scipy import absolute
from scipy.linalg import eig, norm

def pagerank(data):
  n = len(data)
  l, v = eig([[1.0 / len(data[i]) if j in data[i] else 0 for i in range(1, n+1)] for j in range(1, n+1)])
  ev = v[:, (l == absolute(l).max()).nonzero()[0][0]]
  return ev / norm(ev, 1)

data = {
  1: [2, 3, 4, 5, 7],
  2: [1],
  3: [1, 2],
  4: [2, 3, 5],
  5: [1, 3, 4, 6],
  6: [1, 5],
  7: [5],
}

for x in pagerank(data): print float(x)

#2349の7行目は何をやっているのかよくわかりません…
__getslice__がオーバーライドされてるんでしょうか…。

とりあえず僕のコードを貼ってみます。
数学的にはさほど難しいことはしていなくて、
単に固有値を求めて絶対値を取ってから
合計でそれぞれを割って合計が1になるように正規化しています。

data = {
    1: [2, 3, 4, 5, 7],
    2: [1],
    3: [1, 2],
    4: [2, 3, 5],
    5: [1, 3, 4, 6],
    6: [1, 5],
    7: [5],
}

m = [[1.0 / len(data[j]) if i in data[j] else 0
      for j in range(1, 8)]
      for i in range(1, 8)]

from numpy import array
from numpy.linalg import eig

(v, d) = eig(array(m))

pagerank = [abs(row[0]) for row in d]
k = sum(pagerank)
pagerank = [x / k for x in pagerank]

print pagerank
# [0.303514376997, 0.166134185304, 0.140575079872,
#  0.105431309904, 0.178913738019, 0.0447284345048,
#  0.0607028753994]

ScalaがJAMA(Java Matrix Package)を使っていて
Javaが自前で行列演算を実装しているというのも
なんだか奇妙な現象ですね…笑

ライブラリを使えということなので、dojo(0.4.3)を。やってることは#2335とほぼ同じです。

see: dojo

(function(data){ with(dojo.math.matrix){
  var size = 0, i, j, l;
  for(i in data) size++;
  var m = zeros(size, size), pr = create(1, size, 1 / size);
  for(i in data) for(j = l = data[i].length; j--;)
    m[data[i][j] - 1][i - 1] = 1 / l;
  for(var diff = 1, nx; diff > ALMOST_ZERO; pr = nx){
    nx = multiply(m, pr);
    for(diff = 0, j = size; j--;) diff += Math.abs(nx[j][0] - pr[j][0]);
  }
  (typeof WSH != 'undefined' ? function($){WSH.echo($)} : confirm)(format(pr));
}})({
 1: [2, 3, 4, 5, 7],
 2: [1],
 3: [1, 2],
 4: [2, 3, 5],
 5: [1, 3, 4, 6],
 6: [1, 5],
 7: [5]
});

Ruby/GSL版(> GSL-1.9)
Cygwin版Rubyで確認．

require 'gsl'

data = {
  1 => [2, 3, 4, 5, 7],
  2 => [1],
  3 => [1, 2],
  4 => [2, 3, 5],
  5 => [1, 3, 4, 6],
  6 => [1, 5],
  7 => [5],
}

mat = GSL::Matrix.zeros(data.size)
data.each do |k, v|
  v.each do |j|
    mat[j-1, k-1] = 1.0 / v.size
  end
end

# 非対称行列の固有値計算はGSL-1.9以上が必要らしい
eigvals, eigvecs = mat.eigen_nonsymmv
eigvec = eigvecs.abs.get_col(eigvals.abs.max_index)

pagerank = (eigvec / eigvec.sum).to_a
puts pagerank

JAMA使用版を書きました。

不勉強にてJAMAのことも、そもそもScalaがJVM上の言語であることも知りませんでした…orz

最初 jakarta-commons Mathを覗いてみたんですけど、Matrixはあっても固有ベクトル計算はサポートされてなくて、奥村先生のアルゴリズム辞典のヤコビ法のやつを使ったんですが、うまくお題の答えにならなくて(なにが間違っていたかは以下のコードを書いてわかりました。)、結局あのコードになったんですけどね。

package doukaku;
import Jama.Matrix;

public class PageRank {

	public static void main(String[] argv) {
		int[][] data = {{ 2, 3, 4, 5, 7 }, 
				         { 1 }, 
				         { 1, 2 }, 
				         { 2, 3, 5 },
				         { 1, 3, 4, 6 }, 
				         { 1, 5 }, 
				         { 5 }};

		double[][] a = new double[data.length][data.length];
		for (int i = 0; i < data.length; i++) {
			for (int j = 0; j < data[i].length; j++) {
				a[data[i][j] - 1][i] = 1.0 / data[i].length;
			}
		}
		double[] eigVec = new Matrix(a).eig().getV().transpose().getArray()[0];
		double norm1 = new Matrix(eigVec, eigVec.length).norm1();
		for (int i = 0; i < eigVec.length; i++) System.out.println(eigVec[i] / norm1);
	}
}

v[:, i] でi列目のスライスを得られるようですけど、このあたりは
にしおさんのコードでは20行目でrow[0]と、0で決め打ちになってますが、
オリジナルのOctaveのコードではfind(abs(diag(D)) == max(abs(diag(D))))と、
計算でだしてますよね。
この部分を推測をまじえ、固有値の絶対値の最大値が存在する位置をインデックス値として
一列分取り出す、と解釈したものがあの7行目です。

原理が理解できてないんで、ここが0でいいのか、きちんと計算すべきなのか、
わかりません。

ついでに、Numeric付属のMLabを使ったバージョンも投稿します。

MLabのeig関数とSciPyのeig関数では戻り値の行と列が逆になっているような気が。

import MLab

def pagerank(data):
  a = [[1.0 / len(data[i]) if j in data[i] else 0 for i in range(1, len(data)+1)] for j in range(1, len(data)+1)]
  d, v = MLab.eig(MLab.array(a))
  ev = v[MLab.nonzero(d == MLab.max(MLab.absolute(d)))[0], :]
  return (ev / MLab.sum(ev)).real

data = {
  1: [2, 3, 4, 5, 7],
  2: [1],
  3: [1, 2],
  4: [2, 3, 5],
  5: [1, 3, 4, 6],
  6: [1, 5],
  7: [5],
}

for x in pagerank(data): print x

参考のサイト「Googleの秘密 - PageRank徹底解説」での求め方では「絶対値が最大の固有値に対応する固有ベクトル」を取り出すようになっています。しかし、当サイト自身にもありますが、現実には該当する固有ベクトルは必ず第一列目に来るようです。そんなものなのでしょうか。

私も詳しくないのですが。

詳しくは参考URLに書いてありますが、このお題のようなリンク構造のグラフは強連結ですので、対応する遷移行列は常に絶対値が最大1の固有値の固有ベクトルを持ちます。

ですのでそんなものだ、と理解しています。

see: アルゴリズム工学特論

GSLHaskellというライブラリを使った

see: GSLHaskell

import Data.Complex
import Data.List

import GSL.LinearAlgebra.Algorithms (eigR)
import GSL.LinearAlgebra.Matrix (fromLists, toColumns)
import GSL.LinearAlgebra.Vector (toList)

testdata = [[2,3,4,5,7],[1],[1,2],[2,3,5],[1,3,4,6],[1,5],[5]]

pageRank = normalize      -- 確率ベクトルへの正規化
         . map realPart   -- 実数部分の取り出し
         . toList         -- GSLモジュールのベクトルデータからリストに変換
         . head           -- 最大固有値に対応する固有ベクトルの取り出し
         . toColumns      -- 固有ベクトルのリストに変換
         . snd            -- 固有ベクトル(行列データ形式)
         . eigR           -- 固有値計算
         . fromLists      -- GSLモジュールの行列データに変換
         . transMat       -- 遷移確率行列へ変換
         . list2mat       -- 隣接リストから隣接行列へ変換

list2mat xs = map (fillup [1..genericLength xs]) xs 
  where fillup xs [] = map (const 0) xs
        fillup (x:xs) yys@(y:ys) | x < y     = 0 : fillup xs yys
                                 | otherwise = 1 : fillup xs ys

transMat = transpose . map (normalize . map fromInteger)

normalize xs = map (/s) xs where s = sum xs

{-
*Main> pageRank testdata
[0.3035143769968052,0.16613418530351437,0.14057507987220447
,0.10543130990415339,0.17891373801916938,4.4728434504792317e-2
,6.070287539936104e-2]
-}

僕も詳しくはないですが、数学的に「固有値を求める」と表現したときには、どういう手順で求めるのかまでは言及しておらず、結果がどういう順番で並んでいるかは定義されていません。

一方、数値計算ライブラリはなんらかのアルゴリズムで実装されているわけで、アルゴリズムによっては結果が必ずソートされた状態で出てくるものもあるかとは思います。現状は「どうも固有ベクトルが固有値でソートされた状態で出てくるライブラリが多そうだ」ということでしょう。

つまり最大固有値の固有ベクトルが1行目にあるかどうかは、固有値の性質ではなくライブラリの性質なので、ライブラリによっては1行目にないかもしれません。

>固有値の絶対値の最大値が存在する位置をインデックス値として
>一列分取り出す、と解釈したものがあの7行目です。

なるほどなるほど。
「最大固有値は(少なくともNumPyを使って固有値を求めた場合は)0番目にある」
と解釈したのが僕のコードです。

最大固有値が何番目にあるのかをチェックするバージョンも作ってみました。

# あと余談ですがnumpyはLAPACKを使っているということがわかったのでタグを付けときます

see: LAPACK

data = {
    1: [2, 3, 4, 5, 7],
    2: [1],
    3: [1, 2],
    4: [2, 3, 5],
    5: [1, 3, 4, 6],
    6: [1, 5],
    7: [5],
}

m = [[1.0 / len(data[j]) if i in data[j] else 0
      for j in range(1, 8)]
      for i in range(1, 8)]

from numpy import array
from numpy.linalg import eig

(v, d) = eig(array(m))

max_v = list(v).index(max(v))
pagerank = [abs(row[max_v]) for row in d]
k = sum(pagerank)
pagerank = [x / k for x in pagerank]

print pagerank

DouKaku?

PageRankの計算

Posted feedbacks - Flatten