nAG最適化コーナー: 有限差分における導関数の「価値」について

この記事はnAG最適化コーナー シリーズの一部です。

導関数は非線形最適化の分野で重要な役割を果たします。 それは大部分のアルゴリズムが何らかの形で導関数の情報を必要とするからです。

この記事では、微分を計算するいくつかの方法を説明し、よく知られている有限差分近似に詳しく焦点を当てています。

またこの記事全体を通して、目的関数$f(x)$は十分滑らかであり制約なしに最小化されると仮定します。

導関数がなぜ重要なのですか？

典型的な非線形最適化ソルバの動作を思い出してください。このようなソルバはステップ毎に解の推定値を徐々に改善します。 各反復で新しい推定値の位置を決定しますが、その際に利用できるのは目的関数の局所的な情報（現在の関数値と導関数）のみです。

導関数はある点での関数の傾きを表現するので、それを用いて、例えば降下する方向へ、新しい（より良い）点を探索するのは適切と言えます。 さらに、導関数がゼロに近い場合、関数は局所的に水平であり、 局所解に達したことを（他の情報とともに）示します。 導関数の誤りがソルバーを誤解させ、失敗させるかもしれない事は容易に想像できます。

それゆえ、可能な限り、導関数を得るための信頼できる方法を持つことが重要です。

導関数を計算する方法

導関数の評価を手計算で行う事ができるのは関数が解析的に表現できる問題に限られます。 小さな問題以外では、それは面倒で間違いやすいものであり、シミュレーション等のより複雑な関数では不可能です。 このような計算における誤りは非常に一般的な事であり、 それを防ぐために「導関数の診断機能」（ソルバーやライブラリルーチンのオプション等で指定できる場合が多い）の利用を推奨します。

有限差分$FD$（「bumping」とも呼ばれる）は、互いに近い2つの点の関数値の変化で関数の傾きを近似します。 導関数が手作業で導き出されていない場合に最も一般的なアプローチです。これについては、後に詳細を説明します。

別の方法として、自動微分（AD）を介して導関数を計算することもできます。関数の計算がどれほど複雑であっても、非常に多くの場合において、その計算は（非常に）長い一連の基本的な演算と既知の数学関数のみで構成されます。 コンピュータが関数のソースコードにアクセスできる場合、 標準的な微分則（例えば、チェーンルール等）を適用して、導関数の評価を行う事ができます。 しかしながら自動微分は、関数の一部または全部がブラックボックスの場合は適用できません。 また、自動微分はそのまま適用すると、遅すぎる、もしくはメモリが足りない等があり得るので、 自動微分の専門家へのアクセス は常に有用です。

導関数を必要とせず且つ内部的にも有限差分を用いて導関数を用いないアルゴリズムも存在します。 このような一連の解法は「DFO（Derivative-free Optimization）」と呼ばれますが、これについては今後の記事で取り上げて行きたいと思います。

有限差分はどの程度困難になり得るのでしょうか？

小さな実験をしましょう： Rosenbrockのバナナ関数 を一次微分を用いる準ニュートン法 （nAGライブラリルーチン名＝E04KYF）を用いて最小化しようとしてみます。

$$f(x_1,x_2) = 100(x_2-x_1^2)^2 + (1 - x_1)^2$$

これを以下の前進有限差分で近似します。

$$\frac{\partial}{\partial x_i} f(x) \approx \frac{f(x+he_i) - f(x)}{h}$$

ここで$e_i$は標準基底ベクトルです。

実験の結果、有限差分の間隔$h$が重要である事が判明しました。

解は $x_1=x_2=1.0$ で、最適関数値は $0.0$ です。 出発点 $[-2.75, 1.5]$ から数回ソルバーを再実行し、関数評価の回数、最終目的関数値、ステータスを示すifailコード、および最適性の尺度として最終地点の真の勾配のノルムを比較します。

$h$	ステータス $ifail$	関数の評価回数	$f(x)$	$\|\nabla f(x)\|$
正確な導関数	最適解 $0$	57	1.2346e-16	4.9014e-07
1e-13	失敗 $10$	3	3.6894e+03	6.7854e+03
1e-12	最適解 $0$	61	8.7011e-20	2.0029e-09
1e-10	最適解 $0$	60	4.6257e-16	6.1138e-07
1e-08	最適解 $0$	63	8.3944e-12	2.5924e-06
1e-07	準最適解 $7$	188	2.5827e-10	1.4908e-05
1e-06	準最適解 $7$	89	9.0222e-08	3.8266e-04
1e-05	準最適解 $3$	88	8.7574e-06	3.5074e-03
1e-04	準最適解 $3$	79	8.1961e-04	3.8085e-02
1e-03	失敗 $10$	3	3.6894e+03	6.7854e+03

上記の結果から、小さ過ぎるかもしくは大きすぎる $h$ の場合、アルゴリズムが失敗することが分かりました。

$h \in [10^{-12}, 10^{-8}]$ の場合、ソルバーは停止基準を完全に満たし、$h \in [10^{-7}, 10^{-4}]$ の場合、解の妥当な近似が得られたが、$h$ が増加すると信頼性が低下します。

個々の評価における勾配計算の精度も注意深く見ると興味深いものがあります。次の表は $h=10^{-6}$ での実行より選択された点を示しており、有限差分近似の絶対および相対精度、正確および有限差分微分のノルム、および両方のベクトル間の角度（度）を示します：

eval	$\| \nabla f(x) - \nabla f_{FD}(x) \|$	$\| \nabla f(x) - \nabla f_{FD}(x) \| / \| \nabla f(x)\|$	$\| \nabla f(x)\|$	$\| \nabla f_{FD}(x)\|$	angle $deg$
1	4.2385e-03	6.2464e-07	6.7855e+03	6.7855e+03	0.0
5	1.8397e-03	1.5818e-05	1.1630e+02	1.1630e+02	0.0
10	6.7621e-04	2.0091e-06	3.3658e+02	3.3658e+02	0.0
20	1.7706e-04	1.3523e-06	1.3093e+02	1.3093e+02	0.0
30	1.0067e-04	1.3220e-05	7.6146e+00	7.6145e+00	0.0
40	1.3910e-04	2.2372e-05	6.2177e+00	6.2176e+00	0.0
50	3.8461e-04	8.3636e-04	4.5986e-01	4.6013e-01	0.0
51	4.0201e-04	6.3257e-04	6.3552e-01	6.3581e-01	0.0
52	4.0894e-04	2.6149e-03	1.5638e-01	1.5607e-01	0.1
53	4.1319e-04	9.2076e-03	4.4875e-02	4.5189e-02	0.3
54	4.1295e-04	1.6156e-01	2.5561e-03	2.8432e-03	6.3
55	4.1304e-04	1.3006e+00	3.1758e-04	1.7052e-04	67.9
56	4.1305e-04	9.9867e-01	4.1360e-04	6.9688e-07	37.8
57	4.1305e-04	9.9998e-01	4.1306e-04	9.0793e-09	39.7
58	4.1305e-04	1.0017e+00	4.1234e-04	9.2678e-07	39.8

微分方程式の絶対精度はアルゴリズム全体でほぼ同じ$$10^{-4}$$ のままですが、アルゴリズムが解に向かって進むにつれて、傾きは減少しています（53回目以降の評価では$10^{-3}$）。 これにより近似の相対的な精度は破壊され、勾配間の角度は急激に拡大します。

アルゴリズムがこのような状況を「進歩無し」と判断し、2度の反復の後にあきらめている事は驚くべきことではありません。

有限差分区間$h$の効果

ここで観察された挙動は理論によって容易に説明する事ができます。 有限差分近似は、テイラーの定理 を根源とします。 テイラーの定理は、点 $x+d$ における十分に滑らかな関数 $f$ は、点 $x$ における微分値および関数値を用いて次のように表すことができると述べている。
$$f(x+d) = f(x) + \nabla f(x)^T d + \frac{1}{2} d^T \nabla^2 f(x) d + \cdots$$

前進方向の有限差分
$$\frac{\partial}{\partial x_i} f(x) = \frac{f(x+he_i) - f(x)}{h} + \delta_h, \qquad \text{where}\quad \delta_h \approx O(h)$$
は高次の項を無視することによって正確な微分を近似することが容易に分かります。 このことにより $h$ 次の打ち切り誤差 $\delta_h$ が発生します。これは、より小さな $h$ が、より良い近似を与えることを示唆しています。

しかしながら、第2の競合要因も考慮する必要があります。関数値 $f(x)$ 及び $f(x+h)$ の評価は有限精度演算で行われるため、マシン精度 $\varepsilon$ の誤差を考慮する必要があります。 したがって、実数演算における桁落ちによる誤差は $O(\varepsilon/h)$ です。

小さ過ぎる $h$ は打切り誤差を小さくしますが、桁落ちによる誤差を増大させます。

最適な間隔

打切り誤差と桁落ちによる誤差双方の適切なバランスが大切です。よくスケーリングされた関数を仮定すると、$\sqrt{\varepsilon}$ の全体近似誤差で、最適な前進方向有限差分間隔 $h_{opt} \approx \sqrt{\varepsilon}$ を推定することができます。 ([1], [2] 参照)

標準の倍精度浮動小数点演算では $\varepsilon = 1.1e-16$ となるので、$h_{opt}$ はおおよそ $1e-8$ です。

この仮定に合わない関数を想像するのは簡単です。例えばそれぞれの変数のスケールが大きく異なる場合があげられますが、このような場合そのような場合には、有限差分間隔 $h$ に対して注意を払う必要があります。

nAG（ライブラリ）は、$h$ の数値的な推定を行うルーチンを提供します。（e04xac 参照）このルーチンは、座標ごとに少なくとも3回の関数評価を必要とし、最良の使用は出発点などの代表点にあり、計算を通して間隔 $h$ を変更せずに維持します。

各座標に対する異なるスケールの効果は、次の二次関数で示されます。
$$f(x_1,x_2) = (x_1-1)^2 + \alpha*(x_2-1)^2$$
ここで $\alpha \in \{1, 100, 10^4, 10^6\}$ です。

以下では、ルーチン e04xac を点 $[10,10]$ で呼び出し、各座標の最適間隔 $h$ を推定しました。 スケールの変化が $h$ の選択にどのように影響するかが示されました。

$\alpha$	1	100	1e+04	1e+06
$h_{opt}$ for $x_1$	1.1941e-06	8.4602e-06	8.4181e-05	8.4177e-04
$h_{opt}$ for $x_2$	1.1941e-06	8.4602e-07	8.4178e-07	8.4176e-07

発生し得る問題

最適な有限差分間隔 $h$ を選択するための努力にもかかわらず、有限差分近似を損なう要因が存在する場合があります。 前進有限差分は、各点の座標ごとに１回の余分な関数評価しか必要としないという明確な利点を持っています。 しかし、実証されているように、アルゴリズムが解に近づくにつれて近似の相対誤差が大きくなり（ $\|\nabla f(x)\|$ が減少する）、収束が完全に停止する可能性があります。 このような場合は中央有限差分
$$\frac{\partial}{\partial x_i} f(x) = \frac{f(x+he_i) - f(x-he_i)}{2h} + \delta_c, \qquad \text{where}\quad \delta_c \approx O(h^2)$$
は、その打ち切り誤差が $h^2$ であるため、より適切かもしれません。またこの際、最適な有限差分間隔を見つけるために、同様の分析が可能です。 適切にスケーリングされた $f$ を仮定すると、$h_{opt}=\varepsilon^{1/3}$ であり、全体近似誤差は $\varepsilon^{2/3}$ （二重計算ではそれぞれ5e-6および2e-11）となります。

中央差分の欠点は、２回の余分な関数評価が必要となる事です。したがって、前進有限差分近似が不十分である場合に、最後の２回の反復についてのみの代替手段となり得るものです。

さらに高次の微分近似もありますが、通常、標準的な最適化の実行に高い負荷を要求するので、ここで詳細は説明しません。

最適な間隔の推定値 $h_{opt}$ は $f(\cdot)$ の評価の精度に大きく依存します。 関数の評価にノイズがある場合（たとえば、複雑なシミュレーションが目的関数の一部として実行され、完全な $\varepsilon$ 精度に達することができない場合）、有限差分近似の精度は非常に急速に低下し、 唯一の救済策は、DFO（Derivative-free Optimization）のようなノイズにあまり敏感ではない特殊なソルバを選択することです。

覚えておくと役立つ事

ここでは有限差分近似でさえ注意する点がある事を示しました。 具体的には、有限差分間隔 $h$ の選択は、精度にとって重要であり、前進有限差分は解近傍では不十分である可能性があるという点です。 理想的には、ソルバー自身が有限差分近似を制御する事かもしれません。これにより、解近傍で中央差分に切り替えるような事が可能です。

いずれにしても、不正確な関数の評価に有限差分が適用される場合には特別な注意を払うべきであり、 またそのような場合には、DFO（Derivative-free Optimization）のような、より適切なアプローチが存在するかもしれないという事です。

---

参考文献:
[1] Gill P E, Murray W and Wright M H $1981$ Practical Optimization Academic Press
[2] Nocedal J and Wright S J $2006$ Numerical Optimization $2nd Edition$ Springer Series in Operations Research, Springer, New York