ペロン＝フロベニウスの定理でページランク

　最近「ペロン=フロベニウスの定理」というのをたまたま見かけた。これを使うとページランクの計算が高速化されそうだったので、実験してみた。

ペロン=フロベニウスの定理とは？

成分が正である実正方行列には唯一つの最大実固有値が存在し、それに対応する固有ベクトルの各成分は厳密に正である[1]。
らしい。証明はまだ読んでないけど、時間のあるとき読む。

ページランクにどう使うか？

もっとも基本的なページランクのモデルを考えます。
 
　上式のように、あるノードのrankはそのノードにリンクを張っているノードの重み付き和とします。これは行列で表現すると、v = Avとなるようなv（つまりAの固有ベクトル）を求めればよいことになります。ただし、Aはグラフの隣接行列です。

　ここまでは昔から知っていましたが、いくつか疑問点がありました。

• |v| = 10^5程度でも行列Aのサイズが10^10となってしまい、簡単には計算できないのではないか？（そのままだとメモリに乗らないし。）行列が疎であることを利用した効率的な計算方法はないのか？
• 考えているグラフは有向グラフだけど正の実数の固有値は存在するのか？（無向グラフの場合は、隣接行列が対称行列なので実数の固有値が存在するし、トレースが0なので正の固有値が存在するんだけど。）

　ペロン=フロベニウスの定理を見たとき、「これか！」とピーンと来ました。

この定理を利用すると、行列Aの最大固有値は正の実数になります。

また、「任意のベクトルに線型写像を施すとそのベクトルは最大固有値に対応する固有ベクトル方向に近づいていく。」という性質を使うと、最大固有値に対応する固有ベクトルは簡単に計算することができます。さらにさらに、この方法で固有ベクトルを計算する場合は、グラフを隣接行列でもつ必要はなく、隣接リストでOKなので空間計算量もO(E)ですみます。

ソースコード

簡単な例題で正しく計算できていることを確認していますが、バグがあるかもしれません。

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

struct edge {
    int from;
    double cost;
};

typedef vector<vector<edge> > adjList;

const double EPS = 1e-9;

/*
 * perform a linear transform nxt := A * cur.
 * notice that A is represented as an adjacency list.
 */
void transform(const adjList &edgesTo, const vector<double> &cur, vector<double> &nxt) {
    int N = cur.size();
    
    double abs = 0;
    for (int i = 0; i < N; i++) {
        nxt[i] = 0;
        for (int j = 0; j < edgesTo[i].size(); j++) {
            const edge &e = edgesTo[i][j];
            nxt[i] += e.cost * cur[e.from];
        }
        abs += nxt[i] * nxt[i];
    }
    
    // normalize a vector nxt
    abs = sqrt(abs);
    for (int i = 0; i < N; i++)
        nxt[i] /= abs;
}

/*
 * calc the distance between two vectors v and w.
 */
double dist(const vector<double> &v, const vector<double> &w) {
    double d = 0;
    for (int i = 0; i < v.size(); i++)
        d += (v[i]-w[i]) * (v[i]-w[i]);
    d = sqrt(d);
    return d;
}

/*
 * display a vector's elements.
 */
void disp(const vector<double> &v) {
    for (int i = 0; i < v.size(); i++)
        cout << i << ": " << v[i] << endl;
}

int main(int argc, char **argv) {
    int N;
    cin >> N;
    adjList edgesTo = adjList(N);

    for (int i = 0; i < N; i++) {
        int deg;
        cin >> deg;
        for (int j = 0; j < deg; j++) {
            int to;
            cin >> to;
            edgesTo[to].push_back((edge){i, 1./deg});
        }
    }

    vector<double> v[2];
    v[0] = vector<double>(N, 1/sqrt(N));
    v[1] = vector<double>(N);

    for (int i = 0; dist(v[0], v[1]) > EPS; i++) {
        transform(edgesTo, v[i&1], v[(i+1)&1]);
        cerr << i << " th iteration: dist = " << dist(v[0], v[1]) << endl;
    }
    
    // page rank values
    cout << "################################################" << endl;
    cout << "page rank values are below:" << endl;
    cout << "################################################" << endl;
    disp(v[0]);

    return 0;
}

入力のフォーマットは以下のようなかんじ。

3      # ノード数
2 1 2  # v0の出次数は2。v1とv2にリンクしている。
1 0    # v1の出次数は1。v0にリンクしている。
1 1    # v2の出次数は1。v1にリンクしている。

ノード数が10^5, 各ノードの出次数を10^2とした場合、6秒くらいで結果が出た。予想してたより固有ベクトルの収束が速かった。これ使ってSNSのソーシャルグラフを使って実験とかしてみると面白そうな気がする。

references

[1] ペロン＝フロベニウスの定理 - Wikipedia