RとMatlabとOctaveと精度

とある論文をもとにして,Rで計算できるようにコードを書いているんだけど,ちゃんと書けたかどうか自信がないw

ので,恥ずかしながら著者の先生に連絡して,検算してもらえないかお願いしてみたところ,快く引き受けてくださったが,何より元のプログラムがFortranで書かれているようで,現在の環境下では動かないとのこと。で,ありがたいことに,現在お使いのMatlabでコードを書いてみましょう,とお願いしてお返事をいただいたのが八月。ありがたいことに,こちらがお渡ししたサンプルコードにくわえてMatlabのコードも送ってくださったのです。

 

ところが,自分の計算結果と検算する間もなく学会シーズンが始まってしまった。半月ほどあいたが,今やっと時間が取れたので検算してみた。

 

案の定,しょうもないミスがあり(行と列を逆にしてた),それを修正するとだいたい合うんだけど,ちょっとやっぱり違うところが出てくる。

それがまあ奇妙なもので,行列は同じ,固有値は同じなんだけど,固有ベクトルが違うというもの。これは正規化するかしないかとか,計算アルゴリズムの問題だよなあ,どうやって検算しようかと頭をひねった。

ということで,まずMatlabのクローンといわれているOctaveをインストール。これをインストールするために,MacPortsをインストールして,失敗して,探しまわったらdmgファイルが落ちていて何じゃそらってなって,という時間を過ごした後,計算してみた。確かにクローンというだけあって,Matlabのコードはほぼコピペで動く。さて固有値・固有ベクトルというところで,同様の問題が出た。固有ベクトルがちょっと合わない。何だこれ,と思ってみたら,固有値分解するプログラムがMatlab(Octave)は二種類あるんですね。

eig関数とeigs関数で,使うアルゴリズムが違うみたい。後者はスパース行列等にも対応しており,計算パッケージとしてはARPACKをベースにしているとのこと。eigs関数の方を使うと,OctaveでもMatlabと近い(それでも完全一致はしない,小数第一桁目までは合う)数値がでる。

Rは基本,eigen関数だけで,LAPACKがベースになっているようだ。そしてRでARPACKを使うにはigraphパッケージを使う必要があり,そこのARPACK関数を使うとある。ところがこのARPACK関数,正方行列を渡したら固有値分解してくれるんじゃなくて,関数の形でベクトルを預けると,一般化固有値問題を解いて指定した数の固有値を出す・・・そんな使い方をするようだ。たしかに社会ネットワークなど,粗なネットワークならそういう使い方もいいだろうけど,こちらはそういう使い方は嬉しくないわけで。

ということで,ここで手をとめることにした。数値計算科学まで首を突っ込んでいられるほど,時間がないんですよw

まあプログラムの概略はわかったし,だいたい自分が間違えてないことが理解できたからいいかな。あとは俺がMatlabを買うか,Rがもっと発展するかですね!