okkyの日記: perlに慣れるとCがつらい… 8
先だっての話でも一部出てきましたが、「総当りで相関係数を求める話」。
データ入力が " でクォートした Tab Separated Values (TSV)フォーマットで、相関を求めたいエントリの名前が x軸、y軸それぞれ1エントリ1行で書かれたファイルになっています。
で、これらを読み込んで相関を求めるプログラム。今まで perl で書いてありました。一応すでに、使わないデータの undef (メモリ解放/swap IO 防止)、平均値など各データ列単位で事前計算できるものは事前計算しておく、などは対応させました。先だっての話の答のおかげで、絶対計算量も 1/2 近くまで減りました(相関係数 rは r(x,y)と r(y,x)で同じ結果になる。式が対称形をしているから)。アルゴリズムレベルでは、もう打てる手は打った感じです。
これ以上速度を向上させようとすると、Cのようなコンパイラ言語で書いて、SSEとかの機能をコンパイラに使っていただかなくてはいけません(gccなら適切なオプションをつけるとできる)。また、処理をマルチスレッド化して、HTやらコアの台数分ぐらいまでは並列処理しなくてはいけない。つーてもこの場合は HT 分の2台だけだが。
.
しかし。
「データファイルを読み込むコードを書くのが面倒で面倒でたまらない」
一行分のデータを読むためのバッファを用意したりとか(一行のサイズが破廉恥にでかいのでまじめにメモリ管理しなくちゃいけない)。
文字列をパースしたりとか。
値を double に変換したりとか。
こう、計算の第一歩に到達するはるか手前で「うわーーーー、イヤダァ」状態です。
やればできるのは判り切っているし、どうやればいいのかもわかるのですが、その処理が面倒でメモリリークとか結構発生させそうだとか、異常系処理(たまに、データが無いときに "" でエントリが終わっている場合があるのでその対処とか)をperlの替わりにやらなくちゃいけないとか、そういうところまで想像がついてしまって… orz。
しかも作ったプログラムは私が使うだけだしなぁ…(「相関を求めてそこからどうするのか」を理解できる人が回りに少なすぎる。すでに4回ぐらい説明しているのだが、全員お地蔵さん状態になっている…。で挙句に判らない、判らないとうめきながら何日も膨大なデータを目と手で追い回し続け、挙句に間違った結論に到達している…)。どんなに遅いと言っても、すでに人の手よりは何万倍も早いし正確だし…。
.
「ここからここまでは perl で書きます。ここからここまでは C で書きます。ここのこれはこういう配列構造に直して引渡ししてください」
みたいな事を記述すると、自動的にインターフェースだけあわせてくれる処理系ってないかなぁ。今ある「perl から呼び出すライブラリをCで書けますが、C側が perl の処理系を理解して延々ガンバってください」的な奴ではなく…。
いや、perl と同じぐらい文字列からデータを切り出す処理が簡単に書けるなら、そいつとC(C++)の組み合わせでもいいんですけれど…
.
ちなみに、HTのもう一方を持っていかれると、計算している間に pysol ができなくなる(重たくて動かない)のもいやな理由のひとつだ、というのはナイショだ。
そういうモジュールを探す (スコア:1)
Cだと
のであれば、そういう計算をやってくれるモジュールを探す。
Python なら、SciPy。
import scipy
a = scipy.array([data_list])
b = scipy.array([another_data_list])
c = scipy.cov(a, b)
see e.g. [python.org]
Perl でも探せば [google.co.jp]あるかと。
Re:そういうモジュールを探す (スコア:1)
と、書いて、え?これは2変数ではなく、多変数の話なの?っと。
N 変数のデータが M組あるとすれば、
>>> X = scipy.array([[1, 2, 3, 4],
[1, 3, 2, 4],
[4, 3, 2, 1]])
>>> scipy.dot(X, X.T)
array([[30, 29, 20],
[29, 30, 21],
[20, 21, 30]])
みたいな感じ?
無駄な計算やってる?
ライブラリがやってることなので知らない。
Re:そういうモジュールを探す (スコア:1)
相関係数を求める事自体は2変数配列を処理すればいいんですが、「その2変数配列の組み合わせ」がたくさんあるんです。
ようするにm要素のデータの組が n セットある。で、2要素間の相関係数を「全部の組み合わせについて求めよ」というのが問題です。もちろん、自分自身との相関は「1.0か NaN」と決まっているのでそれは省略してOK。でもそのままでも m * (m-1) /2 つまり O(m**2)の組み合わせ全部について相関係数を求めなくてはいけません。
従来は CPAN の Math::NumberCruncher というモジュールを使っていました。これは1セットのデータについて相関係数を出してくれるライブラリだけが提供されています。12時間ぐらい動かした結果によれば、全計算が終わるのにおおざっぱに 60日ぐらいかかる計算になります。これはたまらない。
しかし、相関係数って、「各要素のデータセットごとに事前計算しておける」数字がたくさんあります。それぞれの要素の平均値とか、「各要素と平均値との差分の二乗和の平方根」とか。これらを事前に計算しておくと(それぞれの計算量は O(m)なので) O(m**2) の部分の計算量がかなり減ります。事前計算の段階で「これは、結果は NaN になるな」というものを除去すると O(m**2) の段階では m の数自体を減らせます。で、ここまでやった段階で 1日ぐらいで終わるようになったのですが…。(というここまでが前提)
.
判ると思いますが、この手の計算はベクトル計算ですので、SSE系の命令を駆使すると圧倒的に早くなります。GPGPUを使うともっと早くなるんでしょうが残念ながら、それは持っていないので諦めるとして…。大雑把な予想でも4倍までは軽くて、運がよければ12倍ぐらい早くなります。
しかし。そのためには言語を入れ替えなくてはいけません。データファイルをパースする部分を perl の機能に全面的に依存していたので、ここの部分が面倒で面倒で面倒で…
ベクトル計算で高速化するためにはメモリ上の連続性だの alignment があっているだのが必要で、perl の処理の最後だけちょろっと摩り替えるのではむしろ遅くなります。と言うことは最初から C とかで書かなくてはいけないと言うことで…
とまぁ、そういう話です。ちなみに文字列操作とエラーハンドリングの面倒くささゆえに、すでに諦めの境地に
…駄目ジャン
fjの教祖様
Re:そういうモジュールを探す (スコア:1)
> CPAN の Math::NumberCruncher というモジュール
README [perl.org] を見ても、統計関数とか依存性がないので、別のモジュールを探した方がいいんじゃないかな。
SciPy だと、ぐわんばって計算するところは C とか FORTRAN 由来のコードとかで書かれてるから、インタプリタだから遅いという程度の話は回避できます。
統計だけできればいいのなら、R を入れて、Statistics::R [uwinnipeg.ca] が一番近道かもしれない。
前のコードは積和しか出てなかったので、追加。
>>> X = array([[1, 2, 3, 4],
[1, 3, 2, 4],
[4, 3, 2, 1]])
>>> C = scipy.cov(X)
>>> C
array([[ 1.66666667, 1.33333333, -1.66666667],
[ 1.33333333, 1.66666667, -1.33333333],
[-1.66666667, -1.33333333, 1.66666667]])
Covariance までしか出ないので、そこから先は
>>> for i in range(len(C)):
for j in range(i, len(C)):
print i, j, C[i][j]/scipy.sqrt(C[i][i])/scipy.sqrt(C[j][j])
0 0 1.0
0 1 0.8
0 2 -1.0
1 1 1.0
1 2 -0.8
2 2 1.0
Re:perl で前処理をして、バイナリを書き出す (スコア:1)
C 言語で異なったプロセス間で使い回せる型は結局は 8/16/32/64 bit unsigned/signed int と 32/64 bit float です。構造体は所詮それらの塊。
前処理を perl で終わらせて、C の構造体に対応する物を perl で連続して書き出せば C で fread で一気に読み込めます。perl 自体は知りませんが、バイナリのままで数字などを書き出せますよね?同じシステム上だからエンディアンは同じです。一度きりの使い捨てなら 32bit int と 64bit double に限定すれば、アラインメントも回避できるのでは。
CORBA を使う手もアリかとは思いますが、面倒そうです。
Re:perl で前処理をして、バイナリを書き出す (スコア:1)
おしい。
perl はそもそもがネットワーク越しにシステム監視をするツールを簡便に作るために、awk と sed と csh の文法に C の文法を混ぜて作り上げた言語です。このため「エンディアン」は結構「まじめに」作ってあります。つまり、big endian でも little endian でもどちらでも書けるのですが「デフォルト」は存在しない。
もうひとつおまけに言うと、Cのfloat と double のバイナリフォーマットは機種依存性があります。IEEE754を満たす必要があるのは「演算精度」の問題であって、実はメモリ上のバイナリレベルでの互換性はなくても問題はありません。このため、うかつにメモリを吐かせると、酷いことになります。そして x86の内部doubleは「80bit」あったりするのですよ。強引に64bitに落とせますけどね。
これは魅力的な誘惑ですが、これに負けてはいけません。大抵の場合、これに負けると後で痛い目にあいます。
# そしてそれは、大抵、新しいマシンを買ってもらったり、コンパイラが新しくなったりしたときで、
# 上司などの出資者からすれば急いでいるのでやむなく金を払った、という状態のとき。
# なので移植性の無さゆえに仕事が進んでないことを知るととてつもなく不機嫌になります。
# その次の予算を引きずり出すのにどえらく苦労するぐらいに。
# 金づる^H^H^H偉い人を不機嫌にさせないためにも、この誘惑に負けてはいけません。
もう、ね。80x87と Powerアーキテクチャの浮動小数点フォーマットの違いで何度泣かされた事か…。
と言うわけで、おっしゃる手法はわかるのですが、それは私においては禁じ手なのです。
fjの教祖様
API を使えばよろしいかと Re:perl で前処理をして、バイナリを書き出す (スコア:1)
そうだけど、別の方法として、値を ASCII (CSV形式とか)で書き出すことだってできるけど、Perl の C API [namazu.org] を使えばもっとスマートだよね。
Python にも C/C++ API [python.jp] があって、C の既存ライブラリを使ったり、高速化したいとかいうときは常套手段だからね。
SWIG とかは? (スコア:1)
あ、API というのが
だと言ってるのか。
そいうものでもないと思うが。。。
使ったこと無いのでわかんないけど、SWIG [swig.org]とかはどうなんかな?