Showing posts with label 日記. Show all posts
Showing posts with label 日記. Show all posts

Wednesday, November 10, 2010

乱数発生

 Intel Math Kernel Library (MKL)に入ってある乱数発生ライブラリーVector Statistical Library (VSL)を使うことにする。
 Fortranの組込みであるrandom_numberはpseudo-random number generatorを基礎としているらしいが、具体的な乱数発生方法はわからない(ちゃんと調べてない)。

 一方、VSLには流行のMersenne Twister (これはpseudo-random number)やquasi-random number generatorを簡単に選択することができる。

 pseudo vs. quasiにはたくさんの議論があるよう。なんとなく、後者が積分を計算する上でパフォーマンスが良い(収束が速い)というイメージ。

link
http://demonstrations.wolfram.com/MersenneTwisterAndFriends/

Thursday, November 4, 2010

RCI様様

 Intel MKLはreverse communication interface (RCI)を採用しているので、ユーザーが目的関数やヤコビアンの値を与えるサブルーチンを用意することになる。前者があるのは当たり前だが、後者があるのが面倒、なんて思ったりしたが、実はこのRCIはとてもありがたい。非線形最小二乗法で利用するヤコビアンのルーチン(central difference法)の精度がよろしくなくて、収束しない(プログラムがうまく走るか走らないかはヤコビアンに関する収束条件に大きく依存)なんてことが生じてしまう。もし解析的にヤコビアンを求めることができるなら、それを計算するサブルーチンを用意すればいいのだ。ありがたやRCI。

 自分の場合、ヤコビアンはmaximaを使ってシコシコと。なんとも原始的なアプローチ。コンピュータに詳しいなら、mathematicaとFortranを連携させてこんな無駄な作業をしなくて済むんだろうな。憧れるわー。

(追記)
 Trust Region method(信頼領域法)はglobal convergence(大域的収束)が保証されているという。でも、大域的収束の意味って、単に「少なくとも局所最適解に収束する」ことしか意味してない。もちろん、これは望ましい特徴であるのは間違いない。
 解析的にexactlyな値を与えるJacobian用のサブルーチンを用意しても、Intel MKLの1回のループでは非線形方程式を解いてくれない場合がある。
 ということで、少なくとも2通りのアプローチが考えられる。
 (1) 並列化で遺伝的アルゴリズムを採用
 (2) 収束条件を満たさなければ初期条件を変更するというsequentialな形にする
がある。将来的には別の次元(外側のループ)で並列化することを考えてるから、とりあえず後者の方法を採用してみよう。
 というか、Intel MKLの内容って、RCIと言いながら、この種の反復についてはユーザーが勝手にやってくださいねって感じなんだよな。おかしいよ。

Wednesday, November 3, 2010

makefileを作るべきかな

 コンパイル&リンクするファイルの数が増えると、いちいち入力するのが面倒になる。emacsでも、起動時は必ず入力しないといけない。

 近いうちに、これをみて勉強してみよう。