■フォークト関数の計算（その５）
これまで，超幾何関数，放物柱関数，複素誤差関数を用いたフォークト関数の計算法を紹介してまいりましたが，どのような解法であっても一長一短があり，問題に応じてユーザーが使い分けることが必要になってきます．つまり，いかなる場合でも通用する最良たり得る方法は存在しないと考えられます．そこで，今回のコラムでは各種アルゴリズムの性能比較を示したいと思います．
以下では具体例をあげてこれらのアルゴリズムを比較し，それぞれの特徴や利害得失について述べてみます．例としては，
f(x)=1/π∫(0,∞)exp(-σ^2t^2/2-βt)cos{(γ-x)t}dt
において，σ=1,β=1,c=0の場合のフォークト関数値とその微分・積分値，ヤコビアンdf/dσ,df/dβの値を掲げますが，収束判定基準（打ち切り誤差）を一定のもとで，変数ｘの値を変えたときの収束値を比較してみます．
用いたパラメータとアルゴリズムの相性もかかわってくるし，一回のテストだけで断定的な結論を下すことはできません．テストする問題の多様性が十分な数値実験がなされたのちにだけ，アルゴリズムの優劣を比較することができるのでしょうから，以下の結果は必ずしもアルゴリズムの優劣を示すものではありません．しかし，それでも，収束の安定性（ロバストネス）などアルゴリズムの性質を知る目安ぐらいにはなるであろうと思われます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】フォークト関数値
ｘ 超幾何関数 クンマー変換 放物柱関数 複素誤差関数
0 .208709 .208709 .208709 .208709
1 .165796 .165796 .165796 .165796
2 .090715 .090715 .090715 .090715
3 .043385 .043386 .043386 .043386
4 .022814 .022814 .022814 .022814
5 .013885 .013885 .013885 .013885
6 9.33768D-3 9.33768D-3 9.33768D-3 9.33768D-3
7 6.77862D-3 6.77860D-3 6.77795D-3 6.77860D-3
8 5.03191D-3 5.13668D-3 × 5.22479D-3
9 × 4.03053D-3 × ×
10 × 3.24873D-3 × ×
11 × 2.67525D-3 × ×
12 × 2.24184D-3 × ×
13 × 1.90611D-3 × ×
14 × × × ×
フォークト関数の計算プログラムの試運転したところ，困ったことが発生しました．分布の中央部の計算は正確なのですが，裾の部分の計算が発散あるいは異常値に収束してしまうのです（×印）．
これはベキ級数の収束半径が１ということに起因しています．逆分布関数のように定義域が０〜１の関数であれば問題はないのですが，ベキ級数展開では微分・積分計算が簡単になるかわりに，収束半径の問題がつきまといます．
exp(-bt)を展開して，超幾何関数の和としての求めた結果，関数の裾の部分で計算が発散してしまいましたから，そこで，次にcos{(c-x)t}を展開して，放物柱関数経由で超幾何関数の計算を行いました．足し合わせる順番が変わるので安定に収束することを期待したのですが，結果はexp(-bt)を展開した場合と大差ありませんでした．また，複素誤差関数を用いた場合も，性能に大差はなく，収束安定性，総合計算量とも大体同じでした．
４つのプログラムのなかでもっとも収束半径が大きいのは，超幾何関数にクンマー変換を施したものです．クンマー変換は精度を向上させるために行われるものですが，それによって収束半径も拡大されたことになります．
超幾何関数，放物柱関数，複素誤差関数を用いた方法のいずれにおいても，ベキ級数展開していますから，収束半径は１と考えられます．しかし，それは単位円|z|<1で絶対収束という意味であって，パラメータの値によっては条件収束するはずです．したがって，性能は実際に計算してみるまではまで分からないというところがあるのですが，以上のことにより，収束半径は
クンマー変換》超幾何関数＝複素誤差関数〉放物柱関数
の順と考えられました．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】ヤコビアンの値（df/dσ）
ｘ 超幾何関数 クンマー変換 放物柱関数 複素誤差関数
0 -.099109 -.099109 -.099109 -.099109
1 -.031254 -.031254 -.031254 -.031254
2 .034482 .034482 .034482 .034482
3 .027155 .027155 .027155 .027155
4 .010177 .010177 .010177 .010177
5 4.07235D-3 4.07235D-3 4.07235D-3 4.07235D-3
6 1.80005D-3 1.80006D-3 1.79886D-3 1.80006D-3
7 9.18620D-4 9.19108D-4 7.97332D-4 9.21548D-4
8 2.19729D-4 5.20047D-4 × ×
9 × 3.17046D-4 × ×
10 × 2.04567D-4 × ×
11 × 1.38025D-4 × ×
12 × 9.65616D-5 × ×
13 × 6.96093D-5 × ×
14 × × × ×
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【３】ヤコビアンの値（df/dβ）
ｘ 超幾何関数 クンマー変換 放物柱関数 複素誤差関数
0 -.109601 -.109601 -.109601 -.109600
1 -.060630 -.060630 -.060630 -.060630
2 5.03448D-3 5.03448D-3 5.03448D-3 5.03448D-3
3 .022504 .022504 .022504 .022504
4 .017969 .017969 .017969 .017968
5 .012370 .012370 .012370 .012370
6 8.74331D-3 8.74331D-3 8.74356D-3 8.74356D-3
7 6.46134D-3 6.46113D-3 6.99352D-3 6.46103D-3
8 4.77767D-3 4.95890D-3 × 6.16704D-3
9 × 3.92281D-3 × ×
10 × 3.17952D-3 × ×
11 × 2.62868D-3 × ×
12 × 2.22093D-3 × ×
13 × 1.88279D-3 × ×
14 × × × ×
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【４】微分値（df/dx）
ｘ 超幾何関数 クンマー変換 放物柱関数 複素誤差関数
0 0 0 0 0
1 -.073911 -.073911 -.073911 -.073911
2 -.065116 -.065116 -.065116 -.065116
3 -.031015 -.031015 -.031015 -.031015
4 -.012888 -.012888 -.012888 -.012888
5 -6.06553D-3 -6.06553D-3 -6.06553D-3 -6.06553D-3
6 -3.32018D-3 -3.32018D-3 -3.32018D-3 -3.32018D-3
7 -2.02283D-3 -2.02269D-3 -1.49920D-3 -2.02277D-3
8 -1.29415D-3 -1.32695D-3 × -1.89179D-3
9 × -9.18931D-4 × ×
10 × -6.63281D-4 × ×
11 × -4.94723D-4 × ×
12 × -3.78978D-4 × ×
13 × -2.96814D-4 × ×
14 × × × ×
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【５】積分値（−∞〜ｘ）
ｘ 超幾何関数 クンマー変換 放物柱関数 複素誤差関数
0 .5 .5 .5 .5
1 .693560 .693560 .693560 .693560
2 .820959 .820959 .820959 .820959
3 .885114 .885114 .885114 .885114
4 .916717 .916717 .916717 .916717
5 .934508 .934508 .934508 .934508
6 .945914 .945914 .945914 .945914
7 .953885 .953885 .953895 .953895
8 .959791 .959785 .966230 .959776
9 .982600 .964334 × ×
10 × .967953 × ×
11 × .970901 × ×
12 × .973350 × ×
13 × .975417 × ×
14 × × × ×
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【６】収束半径の問題
超幾何関数，放物柱関数，複素誤差関数などの計算法を用いると，数値微分・数値積分でなしに，Voigt関数に関係する諸値が求まります．ただし，収束半径の問題があり，中央からかなり離れた裾の部分では，計算が発散してしまいます．
収束半径の点では，超幾何関数にクンマー変換を用いたものが最も優れていて，いずれの諸値においても，万能でオールマイティーであるという結果が得られました．また，クンマー変換を用いる方法は，計算しやすさを考えるてみても，適当と考えられます．
この収束半径の問題は，正規分布の尺度母数が小さくなる，すなわち，near Lorentzianでは一層顕著になります．しかし，この欠点は絶望的ではなく，プログラムの運用を工夫すればクリアできると思われます．いずれにせよ，各方法の限界に注意しつつ，長所・短所をしっかり把握して頂きたいものです．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【７】フリップ・フロップの問題
この計算ルーチンを実際に非線形最小２乗法組み込んでcurve fittingを行ったところ，異常収束（best fitだが，非負のはずのparameterが負になってしまう）とか，振動・発散して計算がうまくいきませんでした．
これは，いわゆるフリップ・フロップの問題で，目的関数が単峰性を有するとは限らないため，パラメータの値が一義的に決まらず，計算の収束が不安定になったことが原因と推定されました．
すなわち，Voigt関数をf(x,σ,β)とおくと，この関数はGaussianのscale parameter：σとLorentzianのscale parameter：βを二重に含んでいて，それらが強く相関するためにこのようなことが起こったのだと思われます．
Voigt関数のあてはめ計算自体が，極小値が多数あるフリップ・フロップときていますからかなり厄介な問題です．これはどんな方法についてもあてはまることですが，複雑な多峰性関数であれば，よほどうまく初期値を選ばない限り，最小値に到達していなくても局所的な極小値で止まってしまう可能性があります．収束点が大域的最小点であるという保証はなく，局所的な最小点でしかないこともありますし，鞍点（最大値でも最小値でもない停留点をとる点）のこともあります．つまり，収束値は必ずしも最適値であるとはいえません．
この対策としては，
ａ）なるべく真の最小値に近いと考えられる近似値を初期値として選び，いくつかの初期値について計算を繰り返してやってみる，
ｂ）ひとつの方法の収束値を次の初期値として，ほかのアルゴリズムで再計算する
しかありません．
ところが，非線形関数の大域的な振る舞いは非常に複雑であり，一般に多くの極小点をもちますが，それをしらみつぶしに探索することは不可能に近いのです．なぜならば，人間の直観や勘は１変量から高々３変量までの寡変量の世界ではよく働きますが，多変量の世界では直観は働きにくく，４変量以上の多次元データについてはあまり働かないのが通常であって，しかも次元が大きくなるに従い解の配位は非常に複雑となり，われわれが３次元空間でイメージするものとは大きく異なってくるからです．すなわち，非線形関数はその空間が必ずしもなめらかでない可能性があり，空間がなめらかでないとパラメータの値が少し変わっただけで目的関数値が大きく変化することになり，これは目的関数の応答曲面の安定な評価をする上で好ましくない性質です．
非線形最適化問題の解法を導くにあたって，収束の速い方法を開発することにかなりの努力が注がれ，良好な大域的収束性を持つ方法の開発は，むしろ無視されてきました．実は，大域的な最小値への収束を保証するアルゴリズム（悉皆探索）としては，Metyas型ランダム探索法などの確率的探索法があります．しかし，当然のことながら，悉皆探索では出発点を乱数を利用して十分多く生成するか，乱数を使用しない場合はほとんどすべての格子線上をあたって探し出す必要があり，試行錯誤的に出発点を与える方法に比較して極端に収束速度が遅いという難点があります．確率的探索法と降下法を合成することによって，迅速性と確実性を増した悉皆探索を実行しようとする試みもなされていますが，それとて，いくつかの初期値について計算を繰り返してやってみるという原始的な手法と比べて格段に優れているという保証はなく，上述の困難が近い将来に解決される見込みもたっていません．
悉皆探索法について再検討することは価値があると思われますが，もっともこれにはかなりの時間がかかりそうです．今日，悉皆探索は見はてぬ夢にすぎないのですが，将来とも新しい安定な悉皆探索法を作り出すという考えは痴人の夢でしかないのでしょうか．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【謝辞】
鈴鹿高専の奥井重彦先生より，次のような激励を賜りました．
『宮城県立がんセンター・研究所 佐藤 郁郎先生
先生には，複雑な関数に熱心に取り組んでおられるご様子が拝察できます．理論も実験と同じで，繰り返しアプローチする過程で解決が導かれます．当方の能力では，お尋ねにも，正しく答えられませんが，きっと良い成果が得られるものと期待しております．』
『暑い中，佐藤先生が熱心に取り組まれたVoigt関数の解析的表現が導かれ，外山先生のご研究に成果が期待されるとのことで，大変嬉しい限りです．私も正規分布とローレンツ分布の畳み込みである特殊関数，Voigt関数に始めて出会い，多くの学びをさせていただき，感謝しております．』
こちらこそ，どうもありがとうございます．私のほうこそ，いろいろ勉強させていただきました．これらのプログラムを未発表のまま放置しておくことは取り返しのつかない損失になるだろうと大袈裟にとらえ，ここの発表することにしました．このプログラムが一人でも多くの技術者，研究者，学生に役立つことを願っています．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【参考文献】
１．奥井重彦「特殊関数とその応用」森北出版
２．「新数学公式集」丸善
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
