ようするにm要素のデータの組が n セットある。で、2要素間の相関係数を「全部の組み合わせについて求めよ」というのが問題です。もちろん、自分自身との相関は「1.0か NaN」と決まっているのでそれは省略してOK。でもそのままでも m * (m-1) /2 つまり O(m**2)の組み合わせ全部について相関係数を求めなくてはいけません。
perl はそもそもがネットワーク越しにシステム監視をするツールを簡便に作るために、awk と sed と csh の文法に C の文法を混ぜて作り上げた言語です。このため「エンディアン」は結構「まじめに」作ってあります。つまり、big endian でも little endian でもどちらでも書けるのですが「デフォルト」は存在しない。
そういうモジュールを探す (スコア:1)
Ｃだと
のであれば、そういう計算をやってくれるモジュールを探す。
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)
と、書いて、え？これは２変数ではなく、多変数の話なの？っと。
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]とかはどうなんかな？