ニュートン法
Mathematica の利点のひとつは,導関数を記号的に計算できる点である.Method->"Newton"と指定し,関数が明示的に微分可能であれば,自動的に記号的な導関数の計算が行われるのである.一方で,関数が明示的に近似できる形式になっていなければ,Mathematica は有限差分近似を用いてヘッセ行列を計算する.この際,必要な評価回数を最小にする構造的な情報を用いるか,あるいは変数の数値的値を持つヘッセ行列を返す Mathematica の式を指定することができる.
いくつかのユーティリティ関数を含むパッケージをロードする.
| Out[2]= |  |
次は変数の数値的な値を評価するためだけの関数を定義する.
この関数は変数の数値でのみ評価するよう定義されているので,その導関数は記号的には見付けることができない.
次は,
FindMinimumが勾配とヘッセ行列を計算するのに,有限差分を使わなければならない場合に取るステップを示している.
勾配とヘッセ行列の両方が有限差分を使って計算できた場合,ヘッセ行列の誤差が極めて大きい可能性があるので,別のメソッドを使った方がよい.この場合FindMinimumは極めて正確に最小値を見付けるが,導関数の情報が十分ではないので,それも絶対に確かとはいえない.また,関数と勾配の評価数は記号的導関数を自動的に計算した例題におけるよりもはるかに大きくなる.これは,勾配とヘッセ行列のそれぞれを近似するための余分な評価が必要だからである.
勾配(あるいは自動的に計算できる関数)を与えることが可能であるなら,この方法は一般にはるかにうまくいく.勾配はGradientオプションを使って与えることができる.このオプションには「導関数を指定する」ことのできる方法がいくかある.
これは

と

の数値に対する勾配を返す関数を定義する.
| Out[5]= |  |
| Out[6]= |  |
ヘッセ行列を与えるプログラムがあれば,それも与えるとよい.ヘッセ行列はニュートン法でしか使われないので,
のオプションとして与える.
次で

と

の数値に対するヘッセ行列を返す関数を定義する.
| Out[7]= |  |
| Out[8]= |  |
原則として,ニュートン法は記号微分を評価するか有限差分を使うかして計算されたヘッセ行列を使う.このようにして計算されたヘッセ行列は常に定正値行列であり,メソッドの収束は凸関数に依存する.この条件が破られた位置から探索が始まることもよくあるので,アルゴリズムはこのことを考慮に入れる必要がある.
| Out[9]= |  |
極小値に十分近ければ,ヘッセ行列は実際には負定値行列である.
| Out[10]= |  |
ヘッセ行列が定正値行列ではない場合だけにニュートンのステップ式を適用する場合は,関数値の減少を導かないステップ方向にすることができる.
これは

を解くことで求まる方向についての方向微分を計算する.これは正であるので,この方向に移動すると関数値が局所的に大きくなる.
| Out[11]= |  |
直線探索法の収束では,正定値の二次形式
を使ってその方向を計算することが重要である.そこから導かれる探索過程と収束結果が,十分に降下した方向に依存しているためである.「直線探索法」を参照されたい.Mathematica は
が定正値行列になれるほど大きい要素を持つ対角行列
によってヘッセ行列を修正する.このようなメソッドは修正ニュートン法と呼ばれることがある.密なヘッセ行列でも疎なヘッセ行列でも,
に対する修正は [GMW81]で説明してあるようにコレスキー分解の計算過程のどこかで行われる.この修正は
が定正値行列ではないときにのみ行われる.この分解法を独自に使いたければ,LinearSolveからアクセスすることができる.
次は

を使ってステップを計算する.ここで

はヘッセ行列のコレスキー因子が計算される際に決定される.
| Out[12]= |  |
| Out[13]= |  |
(修正)ニュートン法は,その強力さの他,収束率も主要な側面といえる.探索が極小値に十分近付くと,収束は二次(
)収束といわれる.これは
が極小値の点であるならば,
の定数について以下の式
が成り立つということである.
一般にステップのほとんどは極小値に近付くために使われるので,機械精度では,二次収束により常に重大な相違が生まれるわけではない.しかし,根を非常に高精度で求めたいのであれば,収束が速いという理由から,ニュートン法がたいていの場合最善の選択である.
これはニュートン法を使って非常に高精度の解を計算する.精度は機械精度(初期点の精度)から十万桁の最大作業精度まで順応的に高められる.取ったステップを保存するために,
Reapと
Sowが使われる.関数の評価回数と使われたステップの記録およびその出力にカウンタが使われる.
| Out[14]= |  |
オプションWorkingPrecision->prec が使われると,AccuracyGoalとPrecisionGoalのデフォルトは prec/2になる.このため,この例題では少なくとも5万桁までの最小値を見付けることになる.
これは探索が近付いた最小値の位置での記号解を計算する.
| Out[15]= |  |
これは各ステップの最後の探索点から厳密な最小値までの距離のノルムを計算する.
| Out[16]= |  |
ステップ数よりも関数の評価回数の方が多く必要とされた理由は,Mathematica が初期値の精度から,要求された最高のWorkingPrecisionまで適応的に精度を上げるためである.使用される一連の精度は,点が最小値に収束しつつあるという仮定のもとで,最も高価な最終精度での評価回数をできるだけ少なくするように選ばれる.Mathematica が精度を変更したときに,より高い精度で関数を再評価しなければならないこともある.
Out[17]//TableForm= |
| |  |
一般的にいって,精度は誤差スケール
の約2倍である.ニュートン法では,ステップが計算されるとき,誤差スケールが二次収束に従って実質的に2倍になるので,これが当てはまる.
FindMinimumは常に指定された初期値の精度から始まる.そのため,FindMinimumで適応的な精度制御を使いたくなければ,厳密値,あるいは少なくとも最大のWorkingPrecisionを持つ値から始めることができる.
これは,終了まで10万桁精度のみを使って解を計算する(警告:計算が終了するまでには非常に時間がかかる.)
| Out[18]= |  |
この方法では関数の評価回数は減るが,すべて最高の精度で行われるので,適応精度を使った方が時間がかなり節約できる.例えば,先程のコマンドで適応精度を使わなかったら,機械精度で始めた場合に比べて50倍以上もの時間がかかる.
ニュートン法には,「直線探索法」と「信頼領域法」の両方の刻み幅制御が実装されている.デフォルトは,上記の例にも用いられている直線探索法であるが,これらの多くは信頼領域法のアプローチで行われる.このアプローチでは一般に1ステップつきより多くの数値線形代数計算が必要になる.しかし,ステップがうまく制御されるので,少ない反復回数で収束する可能性がある.
制約条件のない問題のパッケージを使って,狭く曲がった谷を持つ古典的Rosenbrock関数を設定する.
| Out[19]= |  |
| Out[20]= |  |
| Out[21]= |  |
上記2つのプロットの比較から,信頼領域法では探索が谷に沿って行われ,結果として直線探索法よりも少ない関数評価回数で収束しているため,信頼領域法の方がステップをよく制御していることが分かる.
次の表はニュートン法で使えるオプションのまとめである.
| | |
| "Hessian" | Automatic | ヘッセ行列式の計算に使う式 |
| "StepControl" | "LineSearch" | ステップの制御方法. , ,Noneのいずれか |
Method->"Newton"のメソッドオプション