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

 

日記
自粛生活日記07/09

だらだら書くよ。飲んでるからね。     自粛生活も三ヶ月ぐらいになりますかね。 …

日記
insideな毎日

授業が始まった。最初は不安もあったし,学生はどう受け止めるんだろう,と思うことも多かったが,開始から …

日記
いよいよだぞぅ

「気を緩めては行けない」らしいけど、新規の罹患者数も2桁、50人以下の日々も続いて、そろそろ出口が見 …