制約条件のない最適化:Wolfram言語で最適化問題を設定する

導関数の指定

関数FindRootはオプションJacobianを,関数FindMinimumFindMaximumFindFitはオプションGradientを,ニュートン法はメソッドオプションのHessianを持つ.これらの導関数はすべて同一の基本的な構造で指定される.次は導関数の計算方法の指定の方法をまとめたものである.

Automatic関数の記号的導関数を求める.記号的導関数が求まらない場合は有限差分近似を用いる
SymbolicAutomaticと同じ.ただし,有限差分を用いる際には警告メッセージを出す
FiniteDifference導関数の近似に有限差分を用いる
expression変数の局所的数値で与えられた expression を使い,導関数を評価する

gradient,Jacobian,Hessian導関数の計算方法

導関数の基本的な指定とは,その計算方法である.しかし,導関数はすべてオプションも取る.オプションはリスト{method,opts}を使って指定できる.次は導関数のオプションのまとめである.

オプション名
デフォルト値
"EvaluationMonitor"None導関数が評価されるたびに変数の局所値で評価される式.一般に記号的評価を避けるために->ではなく:>で指定される
"Sparse"Automatic導関数のための疎な構造.AutomaticTrueFalse,あるいは非零構造を与えるSparseArray
"DifferenceOrder"1有限差分が導関数の計算に使われる際の差分の精度次数

gradient,Jacobian,Hessian導関数を計算する際のオプション

次はこれらがどのように使われるかの例である.

いくつかのユーティリティ関数を含むパッケージをロードする:
以下で,変数の数値についてのみ評価することのみを意図する関数を定義する:

Method->"Newton"とするだけで,FindMinimumlstolメッセージを出す.それは,導関数についての十分な情報がないために最小値を解くことができなかったからである.

次は,FindMinimumが勾配とヘッセを計算するために有限差分を使わなければならない場合に取るステップを示したものである:

次のように導関数を指定する勾配オプションを使うことができる.

勾配の記号式を使ってf[x,y]の最小値を計算する:

記号的な導関数が常に得られる訳ではない.有限差分からの確度がさらに必要であれば,関数の評価回数を多くするという犠牲を払うことで,差分の精度次数をデフォルトの1よりも大きくすることができる.

勾配の計算に二次精度の有限差分を使ってf[x,y]の最小値を計算する:

ヘッセ行列を近似するために使われる勾配の計算には関数の評価が使われるので,関数の評価回数ははるかに多くなる.(与えられた情報からはこのための記号式は得られないので,ヘッセ行列は有限差分で計算される.)

FindMinimumPlotから与えられる関数,勾配,ヘッセ行列の評価回数についての情報は極めて役に立つ.EvaluationMonitorオプションにより,これが可能になっている.次はそれぞれのタイプの評価の回数を単純に数えるだけの簡単な例である.(評価が行われたときの値を集めるために,プロットはReapSowを使って作成される.)

次はカウンタでステップ数と関数の数,勾配,ヘッセ行列の評価回数の記録を取りながら最小値を計算する:

このような方法は,同様の特徴を持つ問題のクラスでどのようなメソッドおよび/またはメソッドのパラメータを使うと最も成功しやすいかを見極めるのに大変役に立つ.

Wolfram言語が関数の記号的な構造にアクセスできるときは,必要に応じて自動的に関数およびその導関数の構造分析を行い,SparseArrayオブジェクトを使って導関数を表す.それに続く数値線形代数で疎な構造が使えるので,これは全体的な探索効率に大きく影響する.Wolfram言語に構造分析ができないときは,一般に構造は密であると仮定される.しかし,導関数の疎な構造がどのようなものか分かっていれば,"Sparse"メソッドオプションでこれを指定し,導関数の計算(有限差分を使うと,評価回数が大幅に減少できる)とそれに続く線形代数の両方で効率を大幅に上げることができる.これはベクトル値変数を扱うときには特に重要である.この側面をうまく表しているのが,非常に疎な構造を持つ拡張Rosenbrock関数の問題である

次はUnconstrainedProblems`パッケージを使ってFindRootで解ける,1000の変数を持った拡張Rosenbrock関数を記号形式で取り出す:
関数の記号形式を使ってこの問題を解く:

このように単純な形式の関数では,関数のベクトル形式を書くのは簡単である.べクトル形式では,たとえ自動コンパイルを使っても,記号形式のときよりも関数の評価が格段に速くなる.

拡張Rosenbrock関数のベクトル形式を定義する.これは評価効率が大変よい:
問題の構造から初期点をベクトルとして抽出する:
評価にベクトル変数とベクトル関数を使って問題を解く:

この関数の解は,評価は速いが,全体としてはより遅くなる.x_Listパターンのためにヤコビ行列が記号分析に対して不透明なので,ヤコビ行列を有限差分で計算しなければならないからである.これは有限差分が遅いというよりも,ヤコビ行列のすべての列を得るために関数を100回評価しなければならないからである.構造が分かっていれば,ヤコビ行列を得るための評価が2回にできる.この関数ではヤコビ行列の構造は極めて単純である.

以下で,拡張Rosenbrock関数のヤコビ行列について非零構造を持つパターンSparseArrayを定義する.(規則中の値を_と指定することで,出力からも分かるように,SparseArrayPatternタイプのテンプレートと解釈される.)
以下で実際のヤコビ行列の構造を知った上で問題を解く.コストが目覚ましく改善されているのが分かる:

疎な構造が与えられると,評価して疎な構造のテンプレートで指定される位置に対応する値になるような記号式で値を計算することも可能である.値は,SparseArrayで定められている位置に直接対応しなければならない(この順序はArrayRulesを使って見ることができる).指標の矛盾しない順序を知る方法のひとつに,行列を2度転置するというものがある.こうするとSparseArrayの指標が辞書順で並べられる.

以下で非零構造の行列を2度転置して指標を並べ替える:
次は非零構造の行列における指標位置に対応するヤコビ行列における非零要素を返す関数を定義する:
その結果得られる疎な記号ヤコビ行列で問題を解く:

この場合,疎なヤコビ行列を使ってもそれほど速くなる訳ではない.これは,ヤコビ行列があまりに疎なために,ヤコビ行列の有限差分近似が2度の関数評価でしか求められないためと,問題が最小値付近で非常にうまく定義されているので,ヤコビ行列での追加の確度がさほどの差を生まないためである.

変数と初期条件

関数FindMinimumFindMaximumFindRootはすべて同じ形式で変数を指定する.関数FindFitはパラメータの指定に同じ形式を使う.

FindMinimum[f,vars]vars で与えられる変数について,f の極小値を求める
FindMaximum[f,vars]vars で与えられる変数について,f の極大値を求める
FindRoot[f,vars]vars で与えられる変数について,根 を求める
FindRoot[eqns,vars]vars で与えられる変数について,方程式 eqns の根を求める
FindFit[data,expr,pars,vars]dataexpr が最高のフィットを与えるようにするパラメータ pars の値を vars の関数として求める

"Find"関数の変数とパラメータ.

リスト varsFindFitでは pars)は個別の変数指定のリストでなければならない.各変数指定は次の形式でなければならない.

{var,st}変数 var が初期値 st を持つ
{var,st1,st2}変数 var が2つの初期値 st1st2を持つ.2番目の初期条件は主軸法と割線法でしか用いられない
{var,st,rl,ru}変数 var が初期値 st を持つ.探索は var の値が区間[rl,ru]の外に出ると終了される
{var,st1,st2,rl,ru}変数 var が2つの初期値 st1st2を持つ.探索は var の値が区間[rl,ru]の外に出ると終了される

"Find"関数における個々の変数の指定.

vars の指定にはすべて,同じ数の初期値が含まれていなければならない.領域の境界が指定されていないときは,有界ではない,すなわち rl=-ru=とみなされる.

ベクトル値と行列値の変数

変数の最も一般的な使用法は,数を表すことである.しかし,変数入力文法はベクトル,行列,高階のテンソルとして扱われる変数もサポートする.一般に,"Find"コマンドでは,変数は,それに与えられた初期条件と同じ矩形構造をした値を取るものと見なされている.FindFitの場合は例外で,現在スカラー変数しか使えない.

これは行列である:
FindRootを使ってAの固有値と,それに対応する正規化された固有ベクトルを求める:

もちろん,これが固有値を求める最高の方法ではないが,どのようにして初期値から変数の次元が拾われるのかは示している. の初期値は1なので,スカラーであると見なされる.一方 に与えられた初期値は長さ3のベクトルなので,こちらは常に長さ3のベクトルであると見なされる.

変数に複数の初期値を使う場合は,それらの値の次元が一定であること,また初期値の各要素が互いに異なっていることが必要である.

以下では各変数に2つの初期条件を用いて,異なる固有値を求めている:

ベクトル値と行列値が取れる変数の利点のひとつは,大規模問題および/または異なるサイズの問題の自動処理に非常に効率がよい関数が書けることである.

以下でUnconstrainedProblemsパッケージのExtendedRosenbrock問題に相当する目的関数を与える関数を定義する.この関数は,値として2行の行列である を想定している:

の構造が正しくなければ関数の値は無意味なので,この定義はその構造を持つ引数に限られる.例えば,関数を任意のパターンx_について定義し,次いで未定義のシンボルxで評価すると(FindMinimumが行うのはこういうことである),予期せぬ無意味な結果が返される.ベクトル値の変数についての関数を使っている場合は,しばしば定義を制限しなければならないことがある.上記の定義は正しい構造の記号的な値は排除しない点に注意のこと.例えばExtendedRosenbrockObjective[{{x11,x12},{x21,x22}}]はスカラー x11, に対して関数の記号的な表記を与える.

以下では,問題のサイズに対して標準的な値があたえられると,FindMinimumを使って問題を解く.異なるサイズの問題を解くときは,他の何も変えずに の値だけを変更することができる:

Wolfram言語がこの関数の記号導関数を得ることができず,正確ではない有限差分を使わなければならなかったので,解はデフォルトの許容範囲に達していない.

ベクトル値あるいは行列値の変数を使う際の欠点は,現時点ではWolfram言語でこれらの記号的な導関数が計算できないという点である.正しい導関数を与える関数を展開するのがそれほど困難ではないこともある.これに失敗したが,本当にこれ以上の確度が必要である場合は,より高階の有限差分を使うことができる.

次はExtendedRosenbrockObjective関数の勾配を返す関数を定義する.勾配は変数の位置に対応する行列を平坦化することで得られるベクトルである点に注意のこと:
次は勾配の記号的値を使って問題を解く:

ヤコビ行列とヘッセ行列の導関数は疎であることが多い.適切な場合に,これらの導関数の疎な構造を指定することもできる.この指定により,全体的な解の複雑さをかなり軽減することができる.

終了条件

数学的には,滑らかな関数の極小値についての十分な条件は極めて分かりやすい.でヘッセ行列が正定値行列ならば, が極小値となる(ヘッセ行列が半正定値行列であることは必要条件である).根に関する条件はさらに簡単である.しかし,その値がせいぜいある程度の精度までしか分かっていないコンピュータ上で関数 を計算しており,実質的にかなり限られた回数の評価しかできない場合,探索が最小値あるいは根に十分近付いた時点を判断し,有限の許容誤差までの解を計算するためには誤差推定を使う必要がある.ほとんどの場合,このような推定で十分だが,場合によっては関数の未解決な微細動作のために,誤差が生じることがある.

許容範囲は探索が終了するまでに,検索が根や極小値にどのくらい近付こうとするかに影響を及ぼす.通常数値を使って計算された場合に見られることだが,関数自体に何らかの誤差が含まれると仮定すると,使っている数の精度の半分以上の精度で最小値の位置を突き止めるのは不可能なのが一般的である.これは,極小値の二次的性質によるものである.放物線の底付近では最小値から移動する際の高さの変化が非常に緩やかである.このため,関数に何らかの誤差ノイズがあると,そのノイズの平方根にほぼ等しい幅に渡って放物線の実際の上昇が隠されてしまうのが一般的である.これは例題を見るとよく分かる.

いくつかのユーティリティ関数を含むパッケージをロードする:
次のコマンドは,次第に範囲を狭めながら関数の最小値を示す一連のプロットを表示する.機械数で計算された曲線は黒で,100桁精度で計算された実際の曲線は青で示される.

一連のプロットから,機械数のほぼ半分かあるいはそれ以下である次数の変化で,関数の誤差が最小値付近の実際の曲線を隠してしまっているのが明かである.この精度で関数のサンプルを見ただけでは,与えられた点が許容範囲近くの関数の極小値の値であるかどうかを確認する術はない.

導関数の値は,記号的に計算された場合には,信頼性がずっと高くなるが,一般には,導関数の値だけに頼るのは十分ではない.探索では,導関数が一般に許容範囲を満足するほど小さいような関数の極小値を求める必要がある.関数の記号的導関数が計算できず,有限差分や導関数を使わない方法を使った場合,解の確度はより低下する可能性があるので注意が必要である.

根の探索もこれと同じ関数の不正確さの影響を受けることがある.一般にそれほど激しくではないが,誤差推定の中には二次型のメリット関数に基づいているものがある.

このような制約のため,関数Findのデフォルトによる許容範囲はすべて,最終的な作業精度の半分に設定されている.関数にどれほどの誤差があるのかにより,これが到達できる場合もあればできない場合もある.しかしほとんどの場合,これは合理的な目標である.AccuracyGoalオプションとPrecisionGoalオプションを使うと,許容誤差を調整することができる.AccuracyGoal->agPrecisionGoal->pg であれば,許容範囲が と定義される.

FindMinimum を与えると,となるような値 を求めようとする.もちろん,最小値の厳密な位置 は分からないので,は推定される.これは通常過去のステップと導関数の値に基づいて行われる.導関数の条件を最小にするために,追加的な必要条件 が加えられる.FindRootでこれに相当する条件は,根において残差が小さくなるというである.

以下では,が少なくとも12桁精度であるか,の許容範囲内に収まることが分かる.目標精度のを意味しているので,式に影響は与えない.これと同じように目標確度をとすることはできない点に注意する.これは常に残差の大きさに使われるからである:
これは,結果が要求した許容誤差を満足したことを示している:
以下は,関数の最小値を8桁確度まで求めようとする.プロットに見られるように,関数内の誤差のためにFindMinimumは警告メッセージを出す:
最小値における値は基本的に機械イプシロンであるが,位置は次数程度までしか見付からないことを示している:

複数次元では状況がより複雑になる.これは,Freudenstein Roth問題の例におけるように,最小値が比較的狭い谷で見付かった場合等では方向によって誤差にばらつきがあるためである.このような探索の場合はよく探索のパラメータがスケールされ,これが誤差推定に影響を及ぼす.しかし,これにも関わらず,最小値の二次形式が現実的に到達し得る許容範囲に影響を与えるのが一般的である.

根あるいは最小値をデフォルトの許容範囲を超えて求めたければ,最終的な作業精度を上げる必要があるかもしれない.これにはWorkingPrecisionオプションを使う.WorkingPrecision->prec とすると,探索は初期値の精度で始り,探索が収束するにつれて prec になるまで適応的に大きくなる.デフォルトではWorkingPrecision->MachinePrecisionとなっているので,機械数が使われる.これだとずっと速い.精度を上げるとかなり時間がかかるようになるが,関数が適切に定義されていればはるかに正確な結果が得られる.非常に高精度の解の場合はニュートン法を勧める.ニュートン法の二次収束率が,最終的に必要なステップ数を大幅に減らすからである.

WorkingPrecisionオプションの設定値を大きくしても,関数の定義に低い精度の数が使われている場合は意味がないことは記憶しておいた方がよい.一般にWorkingPrecision->prec が効果的になるためには,関数の定義に使われた数は厳密数であるか,少なくとも精度が prec でなければならない.可能であれば,関数に含まれる数の精度は,収束するように,SetPrecisionを使って人為的に prec に上げるとよいが,これは常に可能な訳ではない.いずれにせよ,関数と導関数が数値的に評価されると,内部演算が prec 桁精度で実行できるよう,必要に応じて結果の精度は prec に上がる.それでも,根また最小値の実際の精度・確度およびその位置は,関数の確度によって限定される.このことは,FindFitで精度までしか分からないデータを使うときに,特に注意しなければならないことである.

機械数を使って定義された関数である:
作業精度を高くしても,プロットで示されるのと同じエラーが実際の関数にあるため,最小値はうまく解けない.導関数は,前述のように機械精度での計算と矛盾のないよう指定された:
関数が機械数を持たないときに,20桁精度で計算を行う:

WorkingPrecision->prec と指定して,AccuracyGoalPrecisionGoalオプションを明示的に指定しない場合,Automaticのデフォルトの設定は,AccuracyGoal->prec/2PrecisionGoal->prec/2となる.これにより,前述の通り,一般に現実的に期待できる最小の許容誤差が導かれる.

AccuracyGoalPrecisionGoalオプションの設定を明示的に指定せずに,50桁精度で計算を行う:
最小値の値は,デフォルトの25桁の許容誤差よりも実際にずっとよいことが分かる:

以下の表は,精度および確度に影響するオプションの要約である.

オプション名
デフォルト値
WorkingPrecisionMachinePrecision使用する最終作業精度 prec.精度は prec と初期条件の精度のうち小さい方から始まり,prec まで適応的に増加する
AccuracyGoalAutomatic設定 agtola=10-ag により絶対許容誤差を決定する.Automaticのときag=prec/2である
PrecisionGoalAutomatic設定 pgtolr=10-pg により絶対許容誤差を決定する.Automaticのとき pg=prec/2である

関数"Find"の精度・確度オプション

場合によっては,探索がゆっくりと収束することがある.遅い探索が無限に続くのを防ぐために,Findコマンドには終了までに許される最大反復回数(ステップ)が付いている.これはオプションMaxIterationsで制御でき,デフォルト値はMaxIterations->100である.この条件により探索が終了したときは,コマンドによりcvmitメッセージが表示される.

以下でOptimization`UnconstrainedProblems`パッケージのBrownDennis問題を取り出す.
関数が平方和なので,デフォルトのLevenbergMarquardtメソッドで問題を解く:

LevenbergMarquardtメソッドはこの問題でゆっくりと収束している.それは,残差が最小値付近で非零であり,ヘッセ行列の二次形が必要とされるためである.このメソッドは最終的に400ステップ未満で収束するが,収束の速いメソッドを使った方がよいだろう.

大規模計算では,反復限界に達したときに,別のメソッドとともに返されることのある最終探索点を,探索続行の初期条件として使うことができる.

記号的評価

関数FindMinimumFindMaximumFindRootはすべてHoldAll属性を持っているので,その引数の評価に対して特殊な意味を持つ.まず,変数は第2引数によって決定し,局所化される.次に,関数が記号的に評価され,数値計算に適した形式に処理される.最後に,コマンドの実行中に,関数が異なる数値で繰返し評価される.ここでは,これらのステップを示すリストに説明を加える.

Determine variablesこれは第2引数を処理することで実行される.第2引数が正しい形式(変数と初期値のリスト)ではない場合,正しい形式になるよう評価される
Localize variablesBlockTableの場合と同様に,変数に与えられた割当てが"Find"コマンドのスコープを超えたWolframシステムセッションに影響を与えないよう,また,直前の割当てが値に影響を及ぼさないよう,変数の規則を与える(この段階では,変数は評価されてそれ自身になる)
Evaluate the function局所的に未定義の変数の(記号的)値を使って,第1引数(関数または方程式)を評価する.注意:この点は,バージョン5における変更点なので,以前のバージョンで実行したコードには,多少の調整が必要となる可能性がある.該当する関数で,記号評価により関数が意図されたように維持されない,あるいは異様に遅い場合は,変数の数値に対してのみ評価するよう,関数を定義しなければならない.最も簡単な方法は,f[x_?NumberQ]:=definition のようにPatternTest (?)を使って関数を定義することである
Preprocess the function使用するアルゴリズムの決定 (平方和のときはLevenbergMarquardt法等)を助けるために,関数を解析する.可能であれば,高速の数値評価ができるよう関数を最適化・コンパイルする.FindRootでは,この段階にはまず方程式から関数にすることが含まれる
Compute derivatives可能であれば,必要な記号的導関数を計算する.あるいは,有限差分を使って導関数を計算するのに必要な前処理を行う
Evaluate numerically異なる数値で関数(必要なときは導関数も)を繰返し評価する

"Find"コマンドに対する関数処理のステップ

FindFitHoldAll属性を持たないので,その引数はすべてコマンドが始まる前に評価される.しかし,関数評価の代りに,モデル関数,変数,与えられたデータから最小化する関数を構築するということ以外は,上記の段階をすべて使う.

関数が明示的な式ではないが,値がプログラムの実行によって導かれた場合等,記号評価を避けたい場合もあるだろう.以下の例では,記号評価を避ける方法を示す.

簡単な境界値問題を,狙い撃ち法を使って数値的に解こうと試みる:

関数を数値評価するために,このコマンドは失敗する.Blockの内部で評価するとどうなるか見てみる.

FindRootに与えられた関数を未定義の局所値xpで評価する:

もちろん,これは関数に対して意図されたものではない.xpに依存すらしていない.ここでは,xpに対する数値がないために,NDSolveが失敗している.従って,Firstが第1引数のx[1]を返す規則がないために,ReplaceAll (/.)が失敗する.関数はxpが数値を持たない限り意味をなさないので,そのように定義しなければならない.

以下は,t=-1におけるx'[t]に対する数値関数として,値x[1]を返す関数を定義する:

FindRootの外部に簡単な関数の定義があると,その関数が本当に自分の意図したものであるかどうかを確認するために,個別にテストすることができるという利点がある.

fx1のプロットを作成する:

プロットから,根に対する2つのブラケット値を推測することができるので,「ブレント」法を使って迅速かつ正確に問題を解くことができる.

狙い撃ち問題を解く:

記号評価を避けるためだけに関数を定義しなければならないので,記号評価は面倒なように思える.しかし,記号評価をしなければ,Wolfram言語はそのユニークに統合された数値パワーと記号パワーを利用するのが困難となる.記号評価とは,アルゴリズムの決定,導関数の自動計算,自動最適化と自動コンパイル,構造分析等の記号分析から派生する利点を,コマンドが絶えず利用できるということである.