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がもっと発展するかですね!

 

日記
もう騙せない

蕁麻疹が完治、根治、しないんですよねぇ。 薬が切れてからというもの、たま〜に「怪しい感じだな?」と思 …

日記
誕生日メモ

娘が「牛乳買って。飲まないんだけど」という。 なんで飲まないもんを買うんだ(笑)と返事したら、少し困 …

日記
今年の10大ニュース2019

今年も大晦日になりました。今年もいろんな人に遊んでもらって、感謝しています。例によってぼんやり今年を …