Groovyで行列の固有ベクトルを求める

dev java groovy

ただいま、JGGUG温泉合宿2009に参加中です。
昨日は「どう書く？.org」のGroovyカバレッジ率を上げと、SpringSource Tool SuiteのGroovy対応の評価を兼ねて、こちらのお題をやってました。
どう書く？org「PageRankの計算」 http://ja.doukaku.org/48/
Googleで検索結果の表示順を決定するために使われているので有名なPageRankですが、理論的には「Webサイト間の関連を超巨大疎行列で表現し、最大の固有値に属する固有ベクトルを求める」ことで重要な順にソートできるというシンプルなものなのです。

学生時代には毎日のようにやっていた行列計算ですが、今ではJavaで簡単に行列計算のできるライブラリがいろいろ揃っているので驚きました。

初めに、Apache Commons Mathを試してみました。
http://commons.apache.org/math/
こちらの説明によると、EigenDecompositionを使えば固有値と固有ベクトルが求められるとのこと。
http://commons.apache.org/math/userguide/linear.html
しめしめと思って試してみましたが、残念ながらEigenDecompositionは対称行列(行と列を入れ替えても変化しない行列)の対角化にしか対応しておらず、今回のお題である非対称行列には未対応でした。(涙)

次に、JAMA(Java Matrix Package)を試してみました。
http://math.nist.gov/javanumerics/jama/
こちらは非対称行列の対角化にも対応していますし、何より使い方がシンプルでよいです。

結果、出来上がったコードはこちら。
http://ja.doukaku.org/48/lang/groovy/

import Jama.Matrix

def 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],
]

def dimension = data.size()
def mat = new double[dimension][]

for(int i=0; i<dimension; i++){
    mat[i] = new double[dimension]
    def numLinks = data[i+1].size()
    def invNumLinks = 1/numLinks
    for(int j=0; j<numLinks; j++){
        mat[i][data[i+1][j]-1] = invNumLinks
    }
}

def tmat = new Matrix(mat).transpose()
def eigVec = tmat.eig().getV().transpose().getArray()[0]
def norm1 = new Matrix(eigVec, eigVec.size()).norm1()
eigVec = eigVec.collect{ it/norm1 }
println eigVec

オブジェクトの初期化の手順の関係で、先に行と列が入れ替わった行列を作っておいて、後で転置しています。

大規模な行列計算ではメモリへの格納効率や処理速度などいろいろ最適化しないといけませんが、ちょっとした行列計算ならJAMAでも十分可能ですね。