この文書は Why Not Numerical Recipes? の全訳です。 コメントは 中野武雄 (nakano@webmasters.gr.jp) までお願いします。 まあ NR もその後版を重ねてますし、 この文書の出所・文責もよくわからなくなってますんで、 どこまで信じるかは読者にお任せします (日付も入ってないし)。 以前は "Why not use Numerical Recipes" という、 もうちょっと reference に耐える文書があったように記憶しているのですが。
武井伸光さんには、訳文について多くの有益なコメントをいただきました。
件のページは http://math.jpl.nasa.gov/nr/ ではないか、 というご指摘を関根達夫さんからいただきました。 ただ残念ながらこのページは現在アクセスできないようです (のでリンクにはしてません)。 NR の著者たちからの反論 も見つかりました。いつかこれも訳してみようかと思ってます。
良いニュースは、 このシリーズは理工学で必要となる数値計算のトピックを、 非常に広い範囲に渡ってカバーしていることだ。しかも値段はとても安い。 悪いニュースは、含まれている数学の説明やコードの質や信頼性に、 出来不出来があることだ。 この本に書かれている議論を権威あるものとみなしたり、 このコードの結果の正当性を盲信したりすることは、 危険であると私たちは考えている。
この本の著者たちは、書籍のカバーで「指導的な科学者たち」 と紹介されている。そうでないとする理由は私たちにはない。 しかしながら、彼らが数値解析や数学ソフトウェアに関して、 すばらしい能力を持っていると言うこともできない。 少なくとも、我々がこの本のいくつかの部分を詳細に調べたところでは、 彼らはそのような力量を示しているとは言えない。
この本を扱った、公開されている書評は、2 つのカテゴリに分類できる。 科学者 (ノーベル賞受賞者の Kenneth Wilson も含む) や工学者による推薦や書評においては、 この本の幅広い視野や利便性を、 品質についてはまじめに評価することなく、 激賞する傾向がある。 一方、数値解析の専門家による書評は、 議論やコードの品質に対して非常に批判的だ。
数値解析の専門家によるレビューを 2 つ紹介しよう。
この章における ODE の数値解析法の解説は、視点が 1970 年代のものだ。 もし著者らがこの話題の専門家に助言を求めていたり、 出版されている良質な総説記事に当たっていれば、 彼らはこれらの手法を違った形で評価していただろうし、 これらの手法についても、もっと新しい版を提示していただろうと私は考える。
また Shampine 教授は、数値積分問題における適応的解法が NR で扱われていないことにも言及している。 適応積分は数値解析では非常に評価されている手法であるにもかかわらず、である。
Hanson 博士は、以前 Association for Computing Machinery Transactions on Mathematical Software (ACM TOMS) のアルゴリズム部門の編集者であった。 博士は NR の非線形最小自乗法のコードをテストし、 より知られているコードである MINPAC の LMDIR や NL2SOL の、 公表されている結果と比較を行った。 NR は場合によっては、比較対象のコードに比べ、 20 倍もの繰り返しが必要であることを博士は発見した。 彼のコメントによれば、 Levenberg-Marquardt 法の減衰パラメータλの制御が十分には洗練されていないため、 λのオーバーフローやアンダーフローが起きてしまっているとのことだ。 NR のアルゴリズムは、Marquardt の 1963 年の論文のアイデアを、 骨子に限って実装したものになっている。 このアイデアにはその後 27 年に渡って重要な進歩があった。 我々は LMDIR や NL2SOL、およびそれらの後継コードのほうが、 ずっと効率的で信頼できる [強調編者] と考える。
現在我々が NR の成果に対して注目するようになったのは、 相談を受けたのがきっかけだった。前述の 2 つの話題もその一部である。 そのような相談は他にもあり、その結果我々は 6.6 節の "Spherical Harmonics" と 14.6 節の "Robust Estimation" とを詳細に調べることになった。
6.6 節で与えられている議論・アルゴリズム・コードは、 内部的には整合している。 また Legendre 陪関数の計算に用いている漸化式は、 この話題の専門家たちによって、 安定であると推奨されているものである。 しかしながら、符号や正規化因子に関して、 様々な習慣があることについては何の警告もなされていない。 何も考えずに NR のコードの結果を他のソースと組み合わせると、 おそらく正しくない結果を得ることになってしまう。
ロバスト評価 (robust estimation) の節を読んでみたところ、 我々は「ロバストな直線フィット」を示す図 14.6.1(b) に疑いを持った。 これは「最小自乗フィット」の結果よりもずっと良いように見える。
この図に関する我々の疑問を確かめるため、 我々はこの図を拡大して実験点および「フィットされた」線を方眼紙にトレースし、 その結果に関して次のような実験を行った。
我々はまず最小自乗フィットを計算した。 図 14.6.1(b) において、 どのような「ロバスト」な手法が用いられたのは特定できなかった。 しかし、この領域で NR において与えられているコードは L1 フィッティングだけなので、 我々はこのデータに対する「ロバスト」なフィットの例として、 L1 フィットを計算した。 我々は CL1 サブルーチンを用いて L1 フィットの結果を得た。 これは 1980 年に ACM TOMS のアルゴリズム部門から公開されたもので、 結果に関しては正しいと考えて良いだろう。 我々は同時に NR のコード MEDFIT をこのデータに対して適用した。 得られた結果は、CL1 の結果と小数点以下第三桁まで一致した。
予想通り、最小自乗フィットの結果は 14.6.1(b) の見た目の傾向とあまり違わなかった。 しかし L1 フィットの結果はあまり似ていなかった。 NR の図 14.6.1(b) で "fits" と示された線は、 どんな計算で得られた結果でもなく、 単に著者らが、 「ロバスト」フィッティングに対する支持を補強するために 恣意的に引いたもののように思える。 批判的ではない読者は、おそらく 14.6.1(b) を見て、 これが実際のアルゴリズムの性能である、 と間違って思い込んでしまうだろう。
L1 フィッティング問題の目的関数は、 フィット結果をひとつ以上のデータ点の補間とするようなパラメータ値では、 微分可能ではない。 著者らはこの点に関していくらかは気づいているようだが、 このことが解法アルゴリズムにもたらす結果すべてについてはそうではないらしい。 通常、この問題に対する解法では 2 つ以上の点を補間するのに対し、 著者らのアルゴリズムではアルゴリズムのうちの 試行フィットにおいて、最低 1 つのデータ点しか補間しないことが良くあるようだ。 ここでは MEDFIT/ROFUNC コードが失敗するようなデータセットが 簡単に作れることを示すにとどめておく。
ループを引き起こすようなデータセットのひとつは [x = 1, 2, 3; y = 1, 1, 1] である。コードの別の場所でループを引き起こす別のデータセットとして [x = 2, 3, 4; y = 1, 3, 2] がある。 コードは終了するが、大きく間違った結果をもたらすデータセットは [x = 3, 4, 5, 6, 7; y = 1, 3, 2, 4, 3] である。 理論的な基礎に間違いがあるため、 このコードによって得られた特定の結果を、正しいと信ずることは出来ない。 偶然に正しい結果が得られることも、ないとは言わないが。
"Numerical Recipes" から採った情報やコードに対しては、 その正当性のチェックを別途に行うべきである。
私は第一種の Bessel 関数と修正 Bessel 関数の 0 次と 1 次 (J0, J1, I0, I1) を計算するコードを独自にチェックした。 NR で与えられている近似は、 Cecil Hastings によって 1959 年に出版された、 National Bureau of Standards の "Handbook of Special Functions, Applied Mathematics Series 55" にあるものと同じであった。 この近似は正しいが、非常に正確という訳ではない。 6 桁以上を信じることは出来ない。 よってこのコードを「倍精度」で用いることはできない。 特殊関数の近似に関しては、 この 32 年の間にたくさんの仕事がなされてきているのだが。
我々は NR のアルゴリズムやコードのすべてを一つ一つ調査したわけではないし、 NR の各章のすべての説明をそうしたわけでもない (我々には他にすべき、より生産的な仕事がある)。 しかし (相談を受けて) 任意に採った 4 つの領域において、 その 4 つすべてに問題があったのだから、 残りの部分についても信頼することは難しい。
次に、少しですがコードに奇妙な癖があります。 例えば配列へのアクセスが間違っていたり、 内側のループに不必要な IF や MOD 文を入れていたり。
一方で、連立線型方程式を共役傾斜法で解くことに関する彼らの議論からは、 私はかなりのことを学びました。 彼らの実装は好きではありませんが、でも議論は OK でした。
And:
あなたが "Numerical Recipes" について sci.physics にポストした内容は、 私の経験と一致します。 私は "Numerical Recipes" で与えられている情報は、 ある人を問題に巻き込むものでしかない、と考えています。 なぜなら NR を読んだ人は、 これで何が起こっているかを理解した、と思ってしまうからです。 NR の取り柄は、大抵の話題について参考文献を示しているところです。 何度もひどい目にあえば、まっすぐ参考文献に向かうことを学ぶでしょう。
例: 9.5 節では、多項式の零点を探索するための Laguerre の手法が、どんな出発点からでも収束すると述べています。 しかしながら Ralston と Rabinowitz によれば、 これが成立するのはその多項式の解がすべて実数の場合に限られます。 例えば Laguerre の手法を多項式 f(x)=x^n + 1 に用い、 出発地を 0 とすると、f'(0) = f"(0) = 0 であるために、 厄介な羽目に陥ります。
And:
余談ですが、私はちょうど出版社から、 NR の 18 章と思えるようなプレプリントを受け取りました。 離散ウェーブレット変換に関するものです。 この内容は間違っている、とはっきり言えます。 彼の図にある結果は、彼のルーチンでは再構成できないのです。 なぜかはわかりませんが、とにかく動かないのです。 これらのルーチンを使っている人がもしいらっしゃったら−放り投げることです。 これらのルーチンを修正していたり、 他の離散ウェーブレット変換ルーチンを持っている人がいらっしゃったら、 それについて知りたく思います。宜しく。
And:
私の個人的な警告を。 SVDCMP は常に動作するとは限りません。 ある例では、結果が間違うことがわかりました (幸運なことに確認するのは簡単ですが、普通の人はしないですよね)。 私は NR fortran を c に翻訳し、NR の c コードも試してみました。 両者は同じように間違いました。 IMSL と fortran の Linpack で試し、 Linpack を c に翻訳したものでも試してみました。 この 3 つはすべて正しい結果を示しました。
And:
NR が推薦している乱数生成ルーチン RAN1 と RAN2 は、 まじめなアプリケーションでは決して用いるべきではありません。 RAN1 の最高位ビットを用いて長さ 10,000 の 離散ランダムウォーク (同じ確率で±1) を作ると、 分散は約 1500 で、望ましい値である 10,000 よりも遥かに短いです。
両方とも low-modulus 発生器にかき混ぜバッファを備えたもので、 ある場合に最低位ビットを他の low-modulus 発生器と回転させます。 これらの係数はまじめな仕事には短すぎ、 結果として得られる発生器もとても良いとは言えません。
And:
私が話した人々は皆、この本のそれぞれ違う場所に気に入らないところがありました (私が一番嫌いなのは、simulated annealing と巡回セールスマン問題に関する節です。 この問題についてはずっとましなアプローチが存在します)。
And:
ええ、私も数値解析のずぶの素人で、 NR を究極の言葉として教わり (教授達や同僚達からで、彼らは明らかに NR を使ったことがないのに!)、 なぜ NR の QL 分解のルーチン が動作しないのかを解決するのに 数ヶ月も費やしました。 間違っているのは自分だと思っていたのです。
And:
皆さんの「悲劇集」は、FFT コードに関する私のひどい経験を 裏付けてくれるかも知れません。 もしそのような話があったら教えてください!
And:
私は最近 Numerical Recipes in C (版は不明) のバグに出くわしました。 ksprob() のコードを見てください。 多数回の繰り返しをしても収束に失敗すると、0.0 を返すのです。 しかしこの関数をちょっと見れば、 小さな lambda では和は収束しない (lambda は実際小さくないといけないのですが) ことがわかります。 よって正しい返り値は 0.0 ではなく 1.0 であるべきです。
また、この関数は真の (明らかにまだ誰も解いたことがない) 関数に対する漸近的な近似に過ぎず、この手法は小さな lambda においては 恐ろしく不安定であることを、ユーザは警告されてしかるべきだと思います。
我々の、そして他の多くの人々の経験によれば、 数値計算ソフトウェアを信頼できる情報源から取得するのがもっとも良い。 もっとも手軽かつ安価なのは NETLIB で、 ACM Transactions on Mathematical Software のアルゴリズム (すべて査読を経ている) や、 他にも同業者による精密な調査 (ただし論文誌のレフェリー過程とは少々異なる) に耐えた大量のコードを集めたものである。
MEDFIT の問題は、コードよりも理論のほうが重大だ。 コードの一部を指摘して、「そこだ!」と言うのは難しい。 MEDFIT の作者は、勾配 0 の点で残差が最小になると仮定している。 しかし L1 近似では、これは正しくない。 近似関数は C0 の連続性しか持たないから、 勾配が存在しない点が存在しうる (例えばデータ点はすべてそうなる)。 用いるべき基準は、 "The Approximation of Functions -- Vol 1: Linear Theory" by John R. Rice, Addison-Wesley, 1964 (102ff) の 4-4 節にはっきり書かれている。 最小となるのは、任意の摂動がすべて残差を増加させる場合である。 連続かつ連続的に微分可能な形で接近する場合は、 これは勾配が 0 であることと等しいと言って良い。 しかしこのことは L1 フィッティングの接近では成立しない。 よって MEDFIT のまず最初の欠陥は、 近似関数の勾配が 0 で最小値となる、という仮定にある。
二つ目の欠陥は、問題のある理論をアルゴリズムに転換する際に生じている。 この問題を理解するために、 まずパラメータ 1 つの L1 フィッティング問題を考えてみよう。 この場合、期待すべき結果は集合の中央値である。 しかし集合の大きさが偶数の場合は、 中央の 2 値の間の値なら、どれでも等しく良い解である。 しかしながら NR のアルゴリズムでは、 1 つまたは 2 つのデータ点を補間する解しか考慮していない (一般に N パラメータへの L1 フィッティングは N 個のデータ点を補間する)。 よって結果として、MEDFIT は必ず、微分が存在しない点で微分を計算することになる。 実は MEDFIT が実際に必要とするのは符号だけで、 微分が存在しなくても、正しい値が得られる確率は 1/2 である。 しかしデータがもうちょっとあって、MEDFIT が解の近くにきた場合は、 正と負の勾配が交互に与えられることになり、 2 点の間をギッタンバッコン行き来することとなる。 (おそらくはその 2 点のうちの片方が解だろうが、 MEDFIT にはどちらで終了すればいいのかがわからない)。
3 つめの欠陥は、問題のあるアルゴリズムをコードに転換する際に生じている。 これらの (コードが実際のところ何をしたいのかはわからなくても明らかな) 欠陥については既に他の指摘がある。
ところで、Fortran 版と C 版は、しばしば異なった動作をすることがある。 これは Fortran の SIGN 関数の代わりに用いられている C 版での評価式が、 Fortran-77 標準で定められている SIGN 関数のものとは違うからである。