つれづれ思うこと Perl 編 その３

（その２）の続きとして、次は NumPy について。

例えば、筆者の NMR データの解析のためのプログラムには、exp(i H t) という計算が一杯入ってくる。ここで i は虚数、t は時間であるが、H はハミルトニアンと呼ばれる演算子であり、中身は行列である。つまり、e の何とか乗での冪数（べきすう）にあたる箇所に行列が来る。その行列も１万行 x １万列などという大きさの場合がある（NOESY, saturation transfer のシミュレーションなど）。3 の 2 乗は 9 とすぐに計算できるが、3 の行列乗（しかも内部には虚数）なんて、どのようにして計算するのだろう？と思われるかもしれないが、実はこの行列 (i H t) を対角化さえできれば簡単に計算できる。この exp(i H t) は物理では至るところに出てくるので、物理の計算では行列の対角化（つまり、固有値と固有関数の導出）ができるかどうかが重要な鍵であると言える。

これを Python で計算させると、リターンキーを押し込む前に速攻で答が返ってくる。NumPy と SciPy ライブラリーのおかげである。ここが Python の強いところ。このライブラリー部分は C 言語か Fortran かで書かれているらしいが、表面的にはまるで Python そのものが計算しているように見える。

一方、Perl でこれができるだろうか？一応 PDL（Perl Data Language）と呼ばれるライブラリーがあり、最近ではいろいろな計算ができるらしい。しかし、25 年ほど前は exp(i H t) が Perl で簡単に、しかもバナナジュースを飲み終える短時間内に計算できるかどうかはよく分からなかった。ここが個人的には勝敗を分けたような気がする。1996 年頃、何とかこの計算を Perl で行いたく、力づくでテーラー展開を使った。そうです、対角化という高等な操作を使わなくても、ひたすらテーラー展開をして汗水流せば、exp(i H t) を求めることができるのです。高校でよく x が小さい時 (1+x)^n = 1+nx ですよと習う、あの近似式もテーラー展開から来ている。これを１次でストップせずに、２次、３次と続けていくとますます真の値に近づいていくはずであるが、いったいどこまで続ければよいのだろう？

普通はまあ 10 次ぐらいまで計算して足し込めば十分だろうと思いがちである。ところが、どんどん次数を上げていくと、収束したり発散したりを繰り返したのである（行列の中に虚数が含まれているためだろうか？）。よって、変な次数のところでうっかり止めてしまうと、真値からかなり外れた値に落ちてしまう。というわけで、200 次ぐらいまで計算することにした。この 200 がよいのかどうかは今だに分からない。

ややこしそうに見えるが、要は行列の内積（行と列の要素を一つずつ掛け合わせてから足す）をどんどん繰り返していくだけである。これを計算するには、for 文による繰り返しをひたすら書いていけばよい。そして、NMR のシミュレーションでは FID 部分は 1,000 ポイントぐらいあるので、またこれを 1,000 回ぐらい計算するのである。いずれも for 文を使って。ちょっとややこしいのは、実数と虚数を分けたり混ぜたりしながら計算しないといけない点である（虚数 x 虚数 = 実数, 虚数 x 実数 = 虚数となる具合に）。

このようなプログラムを書いていると何故か汗が出てくる。計算するのはあくまで PC であって、自分が計算するわけではないのであるが、for 文をひとつ追加するごとに、まるで自分が激しく運動しているかのような錯覚に陥り、大変だなあと困って汗ばんでしまうのである。そして何とか書き上げ、いざ実行！だが、答が返ってこない。。。。バグっているのか、それとも PC がまだひたすら計算しているのかは？である。試しに１ポイントだけ計算させてみると、10 分ほど経って一応は正しい答が返ってきた。やった～と叫んだが、たった１ポイントだけではスペクトルをシミュレートするには程遠く、ほとんど役立たずである。それでそのまま放り出してやめてしまった。

月日は流れ 22 年ほど経ったある日、そのプログラムをふと思い出し、こわごわ動かしてみた。普通はそんな昔のプログラムが素直に動くことはない。しかし、この Perl はなにやら計算をし始めた。この Perl の後方互換性は大変すばらしい！私がこれまでに作った Perl の中でうまく動かなったのは一度だけで、それも 30 年以上前に作った初期版だけである。たしか、先頭にスペースのある行を split した時に、配列インデックスが 1 から始まったような気がする（今は０）。それで iMac-27 (core-i5, 32 GB memory）のマシンで 30 分後にプロンプトが返ってきた！恐る恐る gnuplot でグラフを表示させてみると、事前に Mathematica で計算したのと全く同じグラフが出てきたのである。これには大変驚きで、プログラムを前に「なんて可愛いヤツなんだ」と感動してしまった。

あまりもの嬉しさに論文に「Mathematica と Perl で図８をシミュレートした」などと書いたら、ある審査員 reviewer が「こいつはバカか？」と思ったのか、論文を reject すると同時に、NumPy と SciPy のプログラムを送ってきた。それでこれを実行してみると、3 秒ぐらいで同じグラフが画面に出てきたので（gnuplot もなしに）、共著者とともにまたまた驚き、なんて速いプログラムなんだ、大蛇、恐るべし！と、これまた感動してしまった。

それからは、もう勝てないものには勝てないと諦め、素直に Python で、しかし泣く泣く計算プログラムを作るようになった次第である。とは言え、もし Perl にもそのような高速の数学ライブラリが揃っていれば、さらにグラフも簡単に作れれば、間違いなく、蛇ではなくリャマをお友達にしている。