Perl で exp(i H t) などの計算を試してはみたものの、やはり 1995 年頃は遅くて使い物にはならなかった。そこで科学計算にはやはり C 言語を使っていた。その時に役立った本が「Mumerical Recipes in C」（の翻訳版）。ここにはいろいろな計算のためのコードが載っていて、それをそのまま使ってかなりの数値計算をすることができる。特に重宝したのが最適化（フィッティング）の処理である。データに非線形的にフィッティングするための方法では Levenberg–Marquardt 法が有名であるが、コードを打ち込むのが大変そうだったので、何かもっと簡単な方法はないかな？とこの本で探してみると「アメーバ法」という方法が紹介されていた。正式には Nelder–Mead 法と呼ぶらしい。コードも数百行ぐらいで閉じているので、頑張って間違わずに打ち込めば、すぐに動く。この方法のよい点はモデル式の微分が要らない、不連続な箇所があってもよいなどで、まあとにかく制限がほとんど無いので、どんな状況にも適用できる。欠点は時間がかかる、間違えた解にトラップされやすいなど。よって、予想できる解がすでにある時には、できるだけそこからスタートするとよいです。
こんな（失礼ながら）子供遊びのような方法で数学的な最適化ができるのか？と思われ勝ちだが、試してみるとたいへん強力であった。例えば、NMR の残余双極子相互作用値を蛋白質の構造に当てはめる場合、400 残基分のデータをフィットするのに今の PC であれば 0.1 秒とかからない。当時でも 3 秒あれば十分であった。そこで、いろいろな初期値を 10 個ぐらい適当に準備して、それらを全部試してみる。 30 秒ぐらいで 10 個のフィッティングがすべて終わるので、その中から最小の二乗偏差を示す結果だけを選ぶ。エラーバーはどうするのかということになるが、これはモンテカルロ・シミュレーションがよい。100 個ほど模擬データを作っておき、全てにアメーバを振りかける。バナナを食べ終わる 10 分後にはすでに終わっており、その結果をエクセルで統計処理すればエラーバー（標準偏差）が出てくる。
結局つい数年前まで、NMR の緩和、残余双極子相互作用、回転拡散の異方性、拡散係数など、ほとんど全ての処理にアメーバを這わせたが、まったく問題なかった。緩和分散曲線もかなりは問題なし。もちろん Levenberg–Marquardt 法を使えば 0.1 秒ほどで終わってしまうが、今の PC をもってすると、それが延びたところで 1 秒ぐらいなので、もはや大した問題ではなくなっている。今は Perl に書き直したコードも使っている。C の 3-5 倍ぐらいかかるが、上記の NMR データの最適化で 1 分以上かかるような処理はあまりない。
(https://github.com/tikegami-tak/cpmg_hsqc) の中の perl/TI_lib_perl/TI_montecarlo.pm の中の lsquare_simplex というサブルーチンがアメーバ法です。
Fortran 一色という友達がおり「何故？」と尋ねてみると、いろいろ試してみたけれど、やはりこれが最速とのこと。特に数学処理では LAPACK ライブラリなどを使うと（もちろん C でも使えます）そこには何十年もかけて多くの人々が努力してきた成果が埋め込まれており、量子力学の計算ではそれを使わないと計算が終わるまでに一生かかるとのこと。NumPy も中ではこの LAPACK や BLAS といったライブラリーが使われているそうで、まあ道理で速いはずです。そこで、Perl でもこのようなライブラリーをさっと使えるようにならないだろうか？そしたら、アメーバからやっと卒業できるかも？一応、CPAN には Perl と LAPACK などをつなぐラッパーがあるようですが、NumPy のような使い方ではないような気がする。
ちなみに、最近はずるをして Mathematica で Levenberg–Marquardt 法をこっそりと使ったりするのであるが（Perl のアメーバには内緒）、Python を見てみると、変数の範囲を限定して Levenberg–Marquardt 法が使えるようになっていた。さらに CPMG 法では各残基ごとに異なる変数と分子全体で共通の変数があり、Mathematica のようにきっちりとまとまり過ぎた関数ではうまく最適化のためのコードを書けない時がある（ちゃんと勉強したら出来るのかもしれないが、私にはダメだった）。そんな変則的な時、ちょっと原始的な Python は非常に便利だった。また、アメーバ法も今はアルゴリズムがさらに進化しており、もっと高速化しているようです。
0 件のコメント:
コメントを投稿