制約条件のない最適化：非線形方程式の解法

はじめにブレント法
ニュートン法アフィン共変ニュートン法
割線法

はじめに

「極小値」を求めることと非線形方程式の集合を解くことの間には密接な関係がある． 個の未知数を持つ 個の方程式があるとする．解 を求めることは，最小値で残差がゼロのときに の平方和を最小にすることに等しいので，「ガウス・ニュートン」法と特に密接な関係がある．事実，極小値探索のためのガウス・ニュートンステップと，非線形方程式のための「ニュートン」ステップは全く同じである．また，滑らかな関数の場合，極小値探索のための「ニュートン法」は非線形方程式のためのニュートン法と同じである．アルゴリズムの多くの側面が同じであることは驚くべきことではないが，大きな違いも存在する．

この他で最小値探索のアルゴリズムに共通なことに，何らかの「刻み幅制御」が必要なことが挙げられる．一般に刻み幅制御は最小化と同じ方法に基づいている．違いは刻み幅制御の方は，通常滑らかな二乗ノルム であるメリット関数に適用されるという点である．

"Newton"厳密なヤコビ行列あるいは有限差分近似を使って局所的線形モデルに基づいたステップを解く
"Secant"導関数を使わず，過去の 回のステップを用い，ヤコビ行列の割線近似を構築する．各次元に2つの初期条件が必要
"Brent"根の囲い込みを保持する一次元のメソッド．根を囲い込む2つの初期条件が必要
"AffineCovariantNewton"最小数のヤコビアン評価を使うように最適化される

FindRootの基本的なメソッド

ニュートン法

非線形方程式のためのニュートン法は線形近似に基づいている．

このため，ニュートン法のステップは単に と設定することで求めることができる．

ニュートン法は，最小化のための「ニュートン」法と同様に，方程式の根付近で二次収束する．ニュートン法はFindRootのデフォルトのメソッドとして使われている．

ニュートン法は「直線探索法」や「信頼領域法」の刻み幅制御に用いることができる．うまくいくと，通常直線探索制御の方が速いが，信頼領域法によるアプローチの方が強力である．

いくつかのユーティリティ関数を含むパッケージをロードする：
次はFindRoot問題としてのローゼンブロック(Rosenbrock)関数である：
これはデフォルトの直線探索アプローチで非線形系の解を求める．ニュートン法はFindRootのデフォルトのメソッドである：

それぞれの直線探索が直線 に沿って始まっている点に注意のこと．これはこの問題に特有の，ニュートンステップ独特の特性である．

次で，ローゼンブロック関数のヤコビ行列とニュートンステップを記号的に計算する：

このステップが点に追加されると，ステップが直線 に動く理由が簡単に分かる．これはこの問題独特の特徴で，ほとんどの関数では見られないことである．「信頼領域法」はニュートンステップが区間の境界線より内側にないとこの方法を試さないので，信頼領域法による刻み幅制御が使われたときは，このような特徴はあまりはっきりとは現れない．

次は，信頼領域法アプローチを使って非線形系の解を求める．この探索はFindMinimumにおけるローゼンブロック目的関数で「ガウス・ニュートン」法を使った探索とほとんど同一である：

ヤコビ行列の構造が疎である場合，Wolfram言語はヤコビ行列の計算にも，必要な数値線形代数の扱いにもSparseArrayオブジェクトを使う．

非線形方程式を解くことが，陰的メソッドで微分方程式を解くというような，より一般的な数値解法の一部になっているときは，初期値が極めてよいことが多く，完全収束が絶対に必要なわけではないことがある．ニュートンステップの計算で最も高く付く部分は，ヤコビ行列を求めて行列を因数分解する部分であることが多い．しかし，根に十分近ければ，確実に収束率に影響するが，ヤコビ行列を何ステップ分か凍結しておくことも可能である．Wolfram言語でこれを行うのには，ヤコビ行列を更新するまでに何ステップ必要かを指定するメソッドオプションの"UpdateJacobian"を使う．デフォルトは"UpdateJacobian"->1で，この場合，ヤコビ行列は1ステップごとに更新される．

次は，ヤコビ行列が3ステップごとにしか更新されないときに，単純な二乗根を求めるのに必要なステップ数，関数の評価回数，およびヤコビ行列の評価回数を示したものである：
次は，ヤコビ行列が各ステップごとに更新されるときに，単純な二乗根を求めるのに必要なステップ数，関数の評価回数，およびヤコビ行列の評価回数を示したものである：

もちろん，簡単な一次元の根を求めるためにヤコビ行列を更新する場合のコストはそれほど高くない．ここで更新が保存してあるのは例示の目的によるものである．

オプション名
デフォルト値
"UpdateJacobian"1ヤコビ行列を更新するまでに必要なステップ数
"StepControl""LineSearch"刻み幅制御のためのメソッド．"LineSearch"，"TrustRegion"，None（これは推奨しない）のいずれか

FindRootにおけるMethod->"Newton"のメソッドオプション

割線法

導関数が記号的に計算できないときは，「ニュートン」法をヤコビ行列の有限差分近似とともに使う．これは時間と信頼性の両方でコストがかかることがある．その代りに，最小化のときと同じように，導関数なしで使えるように特別に設計されたアルゴリズムを使うこともできる．

一次元では，割線法の考え方は連続する2つの探索点の間の直線の傾きを使って，最新の点における導関数の代りにそのステップを計算するというものである．同様に， 次元では， 個の点における剰余の差を使って，ヤコビ行列の近似に類するものが構築される．これは有限差分に似ているが，できるだけよいヤコビ行列の近似を得るために差分区間を小さくしようとするのではなく，実質的に一次元の割線法のように平均微分を使う．まず， の次元すべてで異なる2つの初期値から 個の点が構築される．ステップが取られるにつれ，最小のメリット関数値を持つ 個の点だけが保存される．稀ではあるがステップが共線であり，割線法によるヤコビ行列の近似が特異値になることもない訳ではない．この場合は，異なる点からアルゴリズムが再スタートされる．

このメソッドは各次元で2つの初期値を必要とする．事実，各次元で2つの初期値が与えられると，一次元を除いて割線法がデフォルトのメソッドになる．一次元では，後で説明するように「ブレント(Brent)」法が選ばれる．

いくつかのユーティリティ関数を含むパッケージをロードする：
次はによるローゼンブロック(Rosenbrock)関数の解を示している：

割線法では，「ニュートン」法より多くの剰余関数の評価が必要になる点に注意のこと．しかし，このメソッドは導関数の情報を直接使わずに比較的狭い谷に沿うことができる．

次は有限差分を使ってヤコビ行列を計算したニュートン法で求めたローゼンブロック問題の解を示している：

しかし，有限差分を使ったニュートン法と比べると，剰余関数の評価回数は同じくらいである．より大きな問題における疎なヤコビ行列の場合，割線法では疎性が利用できないため，有限差分ニュートン法の方が一般に効果的である．

ブレント法

実数値関数の簡単な実数解を求めるときは，関数が負から正に，またはその逆方向に軸を横断するという，この問題の特殊な幾何学性を利用することができる．ブレント法[Br02]は，根が常に囲い込まれるように，関数が正である点と負である点とを常に保持している安全な割線法である．指定された任意のステップで，補間された（割線）ステップといずれ収束することが保証されている2分法のどちらでも選べる．

実関数の根を囲い込む実数の初期条件が2つFindRootに与えられると，ブレント法が使われる．このように，一次元で操作していて根を囲い込む初期条件が決定できるのであれば，そうするのも悪くはない．ブレント法はFindRootで使うことのできるもっとも強力なアルゴリズムだからである．

しかし，たとえ非線形方程式と極小化を解く理論のすべてが滑らかな関数に基づいているとしても，ブレント法は十分に強力なので，不連続関数のゼロ交差もうまく見積れる．

いくつかのユーティリティ関数を含むパッケージをロードする：
次は不連続関数の根の探索で使われたステップと関数評価を示している：

根が極めて近くで囲い込まれると，メソッドはあきらめて警告メッセージを出すが，関数の値(0)を見付けることはできない．この強力さは極めて急な連続関数でも見られる．

次は根付近で急速に変化する関数の根の探索に使われるステップと関数評価を示している：

アフィン共変ニュートン法

ここで示すアフィン共変のニュートン法は，NLEQ_ERRアルゴリズム[Deu06]の実装であり，FindRootのメソッドオプションとして使うことができる．このアルゴリズムの仕組みとその利点を説明するために，アフィン共変のニュートンソルバの簡単な例をいくつか紹介する．これが主に適用されるのは，離散化された非線形偏微分方程式を解く場合である．このような方程式系は大きくなる場合があり，計算速度およびこれらの方程式を解くのに必要なメモリ量の両方におけるコストに注意を払うことが大切である．

まず，方程式の根を求める例題から提示する．非線形方程式系を考える．

ex-2-y=0
y2-x=0

変数と初期値が与えられると，FindRootはこの方程式系の解を求める．

FindRootを使って，ゼロ交差を求める：

このチュートリアルを少し簡単にし，全体を通して一貫性を持たせるために，ベクトル引数を取り，すべての方程式を評価するベクトル値関数 f を設定する．

関数を設定する：
変数名と初期ベクトルをグループ化する：
FindRootを呼出し，変数を値に置き換える：

以下のテキストの比較として，非常に正確であることが知られている基準解があると便利である．

高精度を使って基準解を計算する：

次は，このチュートリアルで提示するタイプの解を求める際にFindRootの呼び出しを簡単にするヘルパー関数である．

FindRootのラッパー関数：
根を求める：

FindRootのデフォルトメソッド等であるニュートン法は，関数 f を評価して計算し，関数 f のヤコビアンを評価する必要がある．ヤコビアンの評価は計算的に高価になる可能性がある．ヤコビアンの評価は，ここで説明するタイプのニュートン法では避けることができないが，ヤコビアンの評価の量を最小限にするアルゴリズムは存在する． FindRootで必要となる評価を調べるために，EvaluationMonitorやStepMonitor等の監視関数を利用することができる．

解を求める段階における関数評価の回数，解を得るのに必要とされるステップ数，ヤコビアン評価の回数を監視することは有益である．

FindRootでのさまざまな評価の回数を調べる：

解を求めるために7ステップと10回の関す評価が必要であった． 関数評価の回数を増加させる試行ステップが試されるが拒否される場合もある．拒否されたステップでは，ステップの数は増えないが，関数評価の回数は増える．解を求める過程を完了させるのに，ヤコビアンの評価が7回必要であったことも分かる．

FindRootで求められたデフォルトの解は，基準解に非常に近い：

FindRootにはいくつかのメソッドオプションがあり，その中のひとつが"AffineCovariantNewton"である．アフィン共変ニュートン法は[Deu06]で詳しく説明されている．

FindRootの"AffineCovariantNewton"メソッドでさまざまな評価の回数を調べる：

アフィン共変ニュートン法では，ステップ数と関数の評価回数は増えたが，ヤコビアンの評価回数は少なくなった．ヤコビアンの評価が計算的に効果である非線形偏微分方程式を解くような場合は，ヤコビアンの評価回数を減らすことでこのような方程式系の解を効率的に求めるという点で大きな差が出る．

Compare to the reference solution:

アフィン共変による解の精度はFindRootのデフォルトメソッドの場合には劣るが，それでもPrecisionGoalのデフォルト設定内である．それを見ていく前に，TrueあるいはFalseに設定できる"AffineCovariantNewton"のメソッドオプション"BroydenUpdates"を見てみよう．ブロイデン更新は，ヤコビアンを再び因数分解することなく，その更新された因子を推定する比較的安価な方法である．アフィン共変ニュートンソルバはいつブロイデン更新を使って効率的に計算するかを見積もる．解を求める段階で，非線形性が増加しブロイデン更新が有用ではなくなると，アルゴリズムはヤコビアンを計算するようもとに戻される．実装されたブロイデン更新アルゴリズムはQNERRアルゴリズム[Deu06]である．

ブロイデン更新なしで系の根を求める：

ブロイデン更新をオフにすると，解を求めるのに必要なステップと関数評価の回数は少なくなるが，ヤコビアンの評価回数は多くなる．

基準解と比較する：

求められた解の精度を増加させるという先ほどの問題に戻ろう．精度は，目標精度を上げることで上がる．

PrecisionGoalを指定して，ニュートン解の差分を減少させる：
基準解と比較する：

アフィン共変ニュートン法は，誤差ノルムが制御されたアルゴリズムであり，AccuracyGoalは監視されない． つまり，このアルゴリズムは誤差だけを監視し，剰余は監視しないのである．誤差はPrecisionGoalに関連付けられており，剰余は目標確度に関連付けられている．

メソッドオプション

アフィン共変ニュートン法は，FindRootオプションの他，多数のオプションを取る．これらのメソッドオプションを説明する．まず，利用可能なオプションのリストを提示する．それから例題とこれらのオプションが有用な場合を説明していく．

以下のメソッドオプションはFindRootのMethod"AffineCovariantNewton"で与えることができる．

オプション名
デフォルト値
"BroydenUpdates"Automaticブロイデン更新
"InitialDampingFactor"Automatic最初のステップにおける減衰係数
"InverseJacobian"Automatic既知の場合，ヤコビアンの逆を指定する
"LinearSolver"Automatic線形のソルバ関数あるいはオプションを指定する
"MinimalDampingFactor"Automatic減衰係数をどれだけ小さくできるか

FindRootのMethod"AffineCovariantNewton"に対するオプション

"BroydenUpdates"

ブロイデン更新は，一旦計算されたヤコビアンを再利用する方法である．前に計算されたヤコビアンを再利用する方が，新しいヤコビアンを再計算するよりも効率的である．欠点は，ブロイデン更新は，新しく計算されたヤコビアンほど正確ではない場合があるということである．アフィン共変ニュートン法アルゴリズムには，安全な場合にはブロイデン更新に変更し，ブロイデン更新が役に立たなくなるとヤコビアンの計算に戻すというメカニズムが実装されている．さらに，ブロイデン更新は無効にすることができる．

ブロイデン更新なしで系の根を求める：
"MinimalDampingFactor"

解法のプロセスのステップ において，方程式 が解かれる．ヤコビアン は逆ヤコビアン隣，新しい増分 が得られる．この増分が前の解に加えられて，更新された解 が得られる． は の間の減衰係数であり，取られる刻み幅を制限する．極端な非線形性がある場合，デフォルトの最小減衰係数 は還元することができる．デフォルト値はである．これが必要となることはめったにない．

"InitialDampingFactor"

非線形性が極端すぎないということが分かっている場合は，より大きい初期減衰係数が選ばれる．デフォルトの初期減衰係数は1/100である．選ばれる値は，"MinimalDampingFactor"より大きくなければならないが，1以下のでなければならない．

大きい"InitialDampingFactor"を指定する：

より大きい初期減衰係数ではヤコビアンの評価を保存することができる．

"InverseJacobian"

ニュートンアルゴリズムには，実は，その役割を果たすために逆ヤコビアンが必要である．逆ヤコビアンはLinearSolveFunctionを生成することで計算される．逆ヤコビアンは，既知であったり計算しやすかったりすることがある．このような場合，逆ヤコビアンはメソッドオプションで"AffineCovariantNewton"に渡すことができる．

評価点 ep において評価された逆ヤコビアンを返す関数を構築する：
問題の関数 に対する逆ヤコビアンを構築する：
指定された逆ヤコビアンで根を求める：

逆ヤコビアンを指定してもヤコビアンの評価回数には影響しない．ヤコビアン評価の回数が0と表示されるのは，逆ヤコビアンがFindRootと相互作用する場合の欠点である．

"LinearSolver"

LinearSolveで使用される線形ソルバやオプションを変更して，逆ヤコビアンを構築することができる．これは"LinearSolver"メソッドオプションを使って行うことができる．

LinearSolveで使用する反復的な"Krylov"メソッドを指定する：

"Krylov"のような反復的メソッドがLinearSolveのメソッドとして指定されると，ヤコビアンの因数分解は可能な場合保存されないが，必要なときに毎回新しく計算される．

線形ソルバは再定義することができる．これによりカスタムの線形ソルバを使うことができる．線形ソルバを再定義することで，ヤコビアンを抽出することもできる．

線形ソルバ再定義してヤコビアンにSowを適用する：
ヤコビアンにReapを適用する：