局所的非線形数値最適化
はじめに
制約条件付きの非線形最適化に対する数値アルゴリズムは,大まかに分けると勾配法と直接探索法とに分けられる.勾配法では,第一導関数(勾配)か第二導関数(ヘッシアン)の乗法が使われる.この例にはSQP(逐次二次計画)法,拡大ラグランジュ法,非線形内点法がある.直接探索法では,導関数情報は使われない.この例としては,Nelder-Mead法,遺伝的アルゴリズム,微分進化法,焼きなまし法がある.直接探索法の方が収束が遅い傾向があるが,関数と制約条件のノイズの存在への耐性は強い.
通常,アルゴリズムは問題の局所的なモデルを構築するだけである.また,そのようなアルゴリズムの多くは目的関数の減少,目的関数と制約条件の組み合わせであるメリット関数の減少を要求し,反復プロセスの収束を確実にする.収束した場合,このアルゴリズムは局所的最適値だけを見付けるため,「局所的最適化アルゴリズム」と呼ばれる.Wolfram言語では,局所的最適化問題はFindMinimumを使って解くことができる.
一方,大域的最適化アルゴリズムは,一般に目的・メリット関数の増加と減少を許可することにより,大域的最適値を見付けようとする.このアルゴリズムは通常計算的により高価である.大域的最適値問題は,Minimizeを使って厳密にとくことも,NMinimizeを使って数値的に解くこともできる.
Minimize[{x - y, -3 x^2 + 2 x y - y^2 ≥ -1}, {x, y}]NMinimize[{x - y, -3 x^2 + 2 x y - y^2 ≥ -1}, {x, y}]FindMinimum[{x - y, -3 x^2 + 2 x y - y^2 ≥ -1}, {x, y}]FindMinimum関数
FindMinimumは,制約条件なしおよび制約条件付きの最適化問題を解く.ここでは制約条件付きの最適化だけを扱う.制約条件なしの最適化でのFindMinimumの詳細は,「制約条件のない最適化」を参照のこと.
FindMinimum[{-(100/(x - 1)^2 + (y - 1)^2 + 1) - (200/(x + 1)^2 + (y + 2)^2 + 1), x^2 + y^2 > 3 }, {x, y}]FindMinimum[{-(100/(x - 1)^2 + (y - 1)^2 + 1) - (200/(x + 1)^2 + (y + 2)^2 + 1), x^2 + y^2 > 3 }, {{x, 2}, y}]前の解の点は実際には局所的最小値である.FindMinimumは局所的最小値だけを見付ける.
ContourPlot[-(100/(x - 1)^2 + (y - 1)^2 + 1) - (200/(x + 1)^2 + (y + 2)^2 + 1), {x, -3, 2}, {y, -3, 2}, RegionFunction -> (#1 ^ 2 + #2 ^ 2 > 3&), Contours -> 10, Epilog -> ({Red, PointSize[.02], Text["global minimum", {-.995, -2.092}], Point[{-.995, -1.992}], Text["local minimum", {0.5304, 1.2191}], Point[{1.2304, 1.2191}]}), ContourLabels -> True]Show[{Plot3D[-(100/(x - 1)^2 + (y - 1)^2 + 1) - (200/(x + 1)^2 + (y + 2)^2 + 1), {x, -4, 3}, {y, -4, 3}, RegionFunction -> (#1 ^ 2 + #2 ^ 2 > 3&), PlotRange -> All], Graphics3D[{Red, PointSize[.02], Text["Global minimum", {-.995, -2.092, -230}], Point[{-.995, -2.092, -207}], Text["Local minimum", {0.5304, 1.2191, -93.4}], Point[{1.23, 1.22, -103.}]}]}]FindMinimumのオプション
FindMinimumは次のオプションを取ることができる.
オプション名
|
デフォルト値
| |
| AccuracyGoal | Automatic | 目標確度 |
| Compiled | Automatic | 関数と制約条件を自動的にコンパイルするかどうか |
| EvaluationMonitor | None | f が評価されるときに常に評価する式 |
| Gradient | Automatic | 勾配関数のリスト{D[f,x],D[f,y],…} |
| MaxIterations | Automatic | 使用する最大反復回数 |
| Method | Automatic | 使用するメソッド |
| PrecisionGoal | Automatic | 目標精度 |
| StepMonitor | None | ステップが取られたときに常に評価する式 |
| WorkingPrecision | MachinePrecision | 内部計算で使用される精度 |
FindMinimum関数のオプション
Methodオプションは,最適化問題を解くのに使用するメソッドを指定する.現在制約条件付き最適化に使用できる唯一のメソッドは内点法である.
FindMinimum[{x^2 + y^2, (x - 1)^2 + 2(y - 1)^2 > 5}, {x, y}, Method -> "InteriorPoint"]MaxIterationsオプションは使用する最大反復回数を指定する.制約条件付きの最適化に機械精度が使われるときは,デフォルトのMaxIterations->500が使われる.
StepMonitorが指定されると,内点法の各反復ステップで1度ずつ評価される.また,EvaluationMonitorが指定されると,同値制約または不同値制約の関数が評価されるたびに,これが評価される.
pts = Reap[sol = FindMinimum[{-(100/(x - 1)^2 + (y - 1)^2 + 1) - (200/(x + 1)^2 + (y + 2)^2 + 1), x^2 + y^2 > 3 }, {{x, 1.5}, {y, -1}}, MaxIterations -> 10, StepMonitor :> (Sow[{x, y}])];]; solContourPlot[-(100/(x - 1)^2 + (y - 1)^2 + 1) - (200/(x + 1)^2 + (y + 2)^2 + 1), {x, -3, 3}, {y, -4, 2}, RegionFunction -> (#1 ^ 2 + #2 ^ 2 > 3&), Contours -> 10, Epilog -> ({Red, PointSize[.01], Line[pts[[2, 1]]], Yellow, Point /@ pts[[2, 1]], Blue, PointSize[.02], Point[pts[[2, 1, 1]]]}), ContourLabels -> True]WorkingPrecision->prec は,FindMinimumの計算をすべて精度 prec で実行するよう指定する.デフォルトは precMachinePrecisionである.prec>MachinePrecisionならば,計算すべてで prec の固定精度が使われる.
AccuracyGoalオプションとPrecisionGoalオプションは,次のように使われる.デフォルトではAccuracyGoal->Automaticであり,prec/3に設定される.デフォルトはPrecisionGoal->Automaticであり,-Infinityに設定される.AccuracyGoal->ga はAccuracyGoal->{-Infinity,ga}と同じである.
AccuracyGoal->{a,ga}およびPrecisionGoal->p であるならば,FindMinimumは,実現可能性とKKT(Karush-Kuhn-Tucker) 条件の満足度と相補条件の組合せである残差を,tol=10-ga 以下にしようとする.また,FindMinimumでは現行と次の反復点 x と x+の差分が,終了前に
を満足することが必要である.
sol = FindMinimum[{-(100 / ((x - 1) ^ 2 + (y - 1) ^ 2 + 1)) - 200 / ((x + 1) ^ 2 + (y + 2) ^ 2 + 1), x ^ 2 + y ^ 2 > 3}, {{x, 3 / 2}, {y, -1}}, WorkingPrecision -> 100, Method -> "InteriorPoint"]solExact = Minimize[{-(100/(x - 1)^2 + (y - 1)^2 + 1) - (200/(x + 1)^2 + (y + 2)^2 + 1), x^2 + y^2 > 3 }, {x, y}];sol[[1]] - solExact[[1]]FindMinimumの例題
大域的最小値を見付ける
大域的最小値が必要な場合,NMinimizeかMinimizeを使わなければならない.しかし,FindMinimumの方がNMinimizeとMinimizeよりも速い傾向にあるので,極小値の少ない比較的大きい問題では特に,FindMinimumを使った方が効率的なことがある.FindMinimumは極小値を求めようとするので,極小値が見つかると終了する.
Plot[Sin[x] + .5 x, {x, -10, 10}]FindMinimum[{Sin[x] + .5 x, -10 <= x <= 10}, {x}]FindMinimum[{Sin[x] + .5 x, -10 <= x <= 10}, {{x, -5}}]FindMinimum[{Sin[x] + .5 x, -10 <= x <= 10 && x < -5}, {x}]SeedRandom[7919];
Table[FindMinimum[{Sin[x] + .5 x, -10 <= x <= 10}, {{x, RandomReal[{-10, 10}]}}], {10}]NMinimize[{Sin[x] + .5 x, -10 <= x <= 10}, {x}, Method -> {"RandomSearch", "PostProcess" -> "InteriorPoint"}]Minimax問題を解く
minimax(minmaxとしても知られる)問題は,いくつかの関数の最大値として定義された関数の最小値を見付けるという問題である.つまり,以下のようなものである.
この問題は,一般の制約条件付き最適化手法を使って解くことができる場合が多いが,問題を滑らかな目的関数を持つ問題に再形式化することによりより簡単に解くことができる.特にminimax問題は次の形式
に変換し,FindMinimumかNMinimizeを使って解くことができる.
(*FindMinMax[{Max[{f1,f2,..}],constraints},vars]*)
SetAttributes[FindMinMax, HoldAll];
FindMinMax[{f_Max, cons_}, vars_, opts___ ? OptionQ] := With[{res = iFindMinMax[{f, cons}, vars, opts]}, res /; ListQ[res]];
iFindMinMax[{ff_Max, cons_}, vars_, opts___ ? OptionQ] := Module[{z, res, f = List @@ ff},
res = FindMinimum[{z, (And @@ cons) && (And @@ Thread[z >= f])}, Append[Flatten[{vars}, 1], z], opts]; If[ListQ[res], {z /. res[[2]], Thread[vars -> (vars /. res[[2]])]}]];FindMinMax[{Max[{x ^ 2, (x - 1) ^ 2}], {}}, {x, y}]FindMinMax[{Max[{Abs[2 x ^ 2 + y ^ 2 - 48 x - 40 y + 304], Abs[-x ^ 2 - 3 y ^ 2], Abs[x + 3 y - 18], Abs[-x - y], Abs[ x + y - 8]}], {}}, {x, y}]ContourPlot[Max[{Abs[2 x ^ 2 + y ^ 2 - 48 x - 40 y + 304], Abs[-x ^ 2 - 3 y ^ 2], Abs[x + 3 y - 18], Abs[-x - y], Abs[ x + y - 8]}], {x, 3, 7}, {y, 1, 3}, Contours -> 40, Epilog -> {Red, PointSize[0.02], Point[{4.93, 2.08}]}]FindMinMax[{Max[{Abs[2 x ^ 2 + y ^ 2 - 48 x - 40 y + 304], Abs[-x ^ 2 - 3 y ^ 2], Abs[x + 3 y - 18], Abs[-x - y], Abs[ x + y - 8]}], {(x - 6) ^ 2 + (y - 1) ^ 2 ≤ 1}}, {x, y}]ContourPlot[Max[{Abs[2 x ^ 2 + y ^ 2 - 48 x - 40 y + 304], Abs[-x ^ 2 - 3 y ^ 2], Abs[x + 3 y - 18], Abs[-x - y], Abs[ x + y - 8]}], {x, 3, 7}, {y, 0, 3}, RegionFunction -> (((#1 - 6) ^ 2 + (#2 - 1) ^ 2 <= 1)&), Contours -> 40, Epilog -> {Red, PointSize[0.02], Point[{5.34, 1.75}]}]多目的最適化:目標計画法
多目的計画法には,制約条件を対象としていくつかの目的関数を最小化することが含まれる.1つの関数を最小化する解は通常他の関数を最小化することはあまりないので,一意の最適解は通常ない.
意思決定者は書く目的に対する目標を想定していることがある.その場合,いわゆる目標計画法を適用することができる.
目標計画法問題をモデル化する方法には多数の変形がある.そのひとつに,優先順位に基づいて目的関数を並べ,重要度の低い目的関数の目標からの偏差を最小化使用とする前に,最も重要な目的関数の目標からの偏差の最小化を求めるというものがある.これは辞書的あるいは先取り目標計画法と呼ばれる.
2つ目の変形は,偏差の重み付き総和を最小化するというものである.特に,次の制約条件付最小化問題が解かれる.
ここで
は実数
の正の部分を表す.重み
は相対的な重要性を反映し,そう敵的なスケールおよび単位を考慮に入れるために偏差を正規化する.重みに対して可能な値は達成すべき目標の逆数である.上の問題は,解きやすい問題に再形式化することができる.
3つ目の変形はチェビシェフ目標計画法というもので,偏差の総和ではなく最大偏差を最小化する.これは異なる目的関数の偏差の均衡を保つ.特に,次の制約条件付き最小化問題を解くとする.
(* GoalProgrammingWeightedAverage[{{{f1,goal1,weight1},...},cons},vars]]*)
GoalProgrammingWeightedAverage[{fg : {{_, _}..}, cons_}, vars_, opts___ ? OptionQ] := With[
{res = Catch[iGoalProgrammingWeightedAverage[{Map[(Append@@#&), Thread[{fg, ConstantArray[1, {Length[fg]}]}]], cons}, vars]]},
res /; ListQ[res]
];
GoalProgrammingWeightedAverage[{fg : {{_, _, _}..}, cons_}, vars_, opts___ ? OptionQ] := With[
{res = Catch[iGoalProgrammingWeightedAverage[{fg, cons}, vars]]},
res /; ListQ[res]
];
iGoalProgrammingWeightedAverage[{fg : {{_, _, _}..}, cons_}, vars_, opts___ ? OptionQ] := Module[
{fs, goals, zs, z, res, ws},
{fs, goals, ws} = Transpose[fg];
If[!VectorQ[ws, (# >= 0&)], Throw[$Failed]];
If[!VectorQ[goals, ((NumericQ[#] && Head[#] =!= Complex)&)], Throw[$Failed]];
zs = Array[z, Length[fs]];
res = FindMinimum[{ws.zs, (And@@Flatten[{cons}, 1]) && (And@@Thread[zs ≥ 0]) && (And@@Thread[zs ≥ fs - goals ])}, Join[Flatten[{vars}, 1], zs], opts];
If[ListQ[res], {fs /. res[[2]], Thread[vars -> (vars /. res[[2]])]}]
];(* syntax GoalAttainment[{{{f1,goal1,weight1},...},cons},vars]]*)
GoalProgrammingChebyshev[{fg : {{_, _}..}, cons_}, vars_, opts___ ? OptionQ] := With[
{res = Catch[iGoalProgrammingChebyshev[{Map[(Append@@#&), Thread[{fg, ConstantArray[1, {Length[fg]}]}]], cons}, vars]]},
res /; ListQ[res]
];
GoalProgrammingChebyshev[{fg : {{_, _, _}..}, cons_}, vars_, opts___ ? OptionQ] := With[
{res = Catch[iGoalProgrammingChebyshev[{fg, cons}, vars]]},
res /; ListQ[res]
];
iGoalProgrammingChebyshev[{fg : {{_, _, _}..}, cons_}, vars_, opts___ ? OptionQ] := Module[
{fs, goals, y, res, ws},
{fs, goals, ws} = Transpose[fg];
If[!VectorQ[ws, (# >= 0&)], Throw[$Failed]];
If[!VectorQ[goals, ((NumericQ[#] && Head[#] =!= Complex)&)], Throw[$Failed]];
res = FindMinimum[{y, (And@@Flatten[{cons}, 1]) && (And@@Thread[y ≥ ws * (fs - goals) ])}, Append[Flatten[{vars}, 1], y], opts];
If[ListQ[res], {fs /. res[[2]], Thread[vars -> (vars /. res[[2]])]}]
];res1 = GoalProgrammingWeightedAverage[{{{x ^ 2 + y ^ 2, 0}, {4(x - 2) ^ 2 + 4(y - 2) ^ 2, 0}}, y - x == -4}, {x, y}]res2 = GoalProgrammingChebyshev[{{{x ^ 2 + y ^ 2, 0}, {4(x - 2) ^ 2 + 4(y - 2) ^ 2, 0}}, y - x == -4}, {x, y}]g1 = ContourPlot[x ^ 2 + y ^ 2, {x, 2, 6}, {y, -1, 2}, ContourShading -> False, ContourStyle -> Blue, ContourLabels -> Automatic];
g2 = ContourPlot[4(x - 2) ^ 2 + 4(y - 2) ^ 2, {x, 2, 6}, {y, -1, 2}, ContourShading -> False, ContourStyle -> Red, ContourLabels -> Automatic];
Show[{g1, g2}, Epilog -> {Line[{{3, -1}, {6, 2}}], PointSize[0.02], Yellow, Point[{x, y} /. res1[[2]]], Green, Point[{x, y} /. res2[[2]]]}]適用例:ポートフォリオ最適化
投資管理のパワフルな方法は,相関関係がほとんど,あるいは全くない資産へ投資することでリスクを分配することである.例えば,ある年に資産Aは20%増加して翌年10%減少し,資産Bはある年に10%減少して翌年20%増加するというふうに,Aの上昇年はBの下降年とすると,両方を同じだけ保有していると,リスクなしで毎年10%の増加になる.実際にはそのような資産はほとんどあり得ないが,この考え方は役に立つものである.
この例の目的は,株,債権,金に投資して,リスクを最低限に抑え,事前設定したリターンを得るような最適資産分配を見付けるというものである.
次は1973年から1994年までのさまざまな資産のリターン歴である.例えば,1973年にはS&P 500は1-0.852=14.8%減少しているが,金は67.7%増加している.
| "3m Tbill" | "long Tbond" | "SP500" | "Wilt.5000" | "Corp. Bond" | "NASDQ" | "EAFE" | "Gold" | |
| 1973 | 1.075` | 0.942` | 0.852` | 0.815` | 0.698` | 1.023` | 0.851` | 1.677` |
| 1974 | 1.084` | 1.02` | 0.735` | 0.716` | 0.662` | 1.002` | 0.768` | 1.722` |
| 1975 | 1.061` | 1.056` | 1.371` | 1.385` | 1.318` | 0.123` | 1.354` | 0.76` |
| 1976 | 1.052` | 1.175` | 1.236` | 1.266` | 1.28` | 1.156` | 1.025` | 0.96` |
| 1977 | 1.055` | 1.002` | 0.926` | 0.974` | 1.093` | 1.03` | 1.181` | 1.2` |
| 1978 | 1.077` | 0.982` | 1.064` | 1.093` | 1.146` | 1.012` | 1.326` | 1.295` |
| 1979 | 1.109` | 0.978` | 1.184` | 1.256` | 1.307` | 1.023` | 1.048` | 2.212` |
| 1980 | 1.127` | 0.947` | 1.323` | 1.337` | 1.367` | 1.031` | 1.226` | 1.296` |
| 1981 | 1.156` | 1.003` | 0.949` | 0.963` | 0.99` | 1.073` | 0.977` | 0.688` |
| 1982 | 1.117` | 1.465` | 1.215` | 1.187` | 1.213` | 1.311` | 0.981` | 1.084` |
| 1983 | 1.092` | 0.985` | 1.224` | 1.235` | 1.217` | 1.08` | 1.237` | 0.872` |
| 1984 | 1.103` | 1.159` | 1.061` | 1.03` | 0.903` | 1.15` | 1.074` | 0.825` |
| 1985 | 1.08` | 1.366` | 1.316` | 1.326` | 1.333` | 1.213` | 1.562` | 1.006` |
| 1986 | 1.063` | 1.309` | 1.186` | 1.161` | 1.086` | 1.156` | 1.694` | 1.216` |
| 1987 | 1.061` | 0.925` | 1.052` | 1.023` | 0.959` | 1.023` | 1.246` | 1.244` |
| 1988 | 1.071` | 1.086` | 1.165` | 1.179` | 1.165` | 1.076` | 1.283` | 0.861` |
| 1989 | 1.087` | 1.212` | 1.316` | 1.292` | 1.204` | 1.142` | 1.105` | 0.977` |
| 1990 | 1.08` | 1.054` | 0.968` | 0.938` | 0.83` | 1.083` | 0.766` | 0.922` |
| 1991 | 1.057` | 1.193` | 1.304` | 1.342` | 1.594` | 1.161` | 1.121` | 0.958` |
| 1992 | 1.036` | 1.079` | 1.076` | 1.09` | 1.174` | 1.076` | 0.878` | 0.926` |
| 1993 | 1.031` | 1.217` | 1.1` | 1.113` | 1.162` | 1.11` | 1.326` | 1.146` |
| 1994 | 1.045` | 0.889` | 1.012` | 0.999` | 0.968` | 0.965` | 1.078` | 0.99` |
| average | 1.078 | 1.093 | 1.120 | 1.124 | 1.121 | 1.046 | 1.141 | 1.130 |
R = {{1.075, 1.084, 1.061, 1.052, 1.055, 1.077, 1.109, 1.127, 1.156, 1.117, 1.092, 1.103, 1.08, 1.063, 1.061, 1.071, 1.087, 1.08, 1.057, 1.036, 1.031, 1.045}, {0.942, 1.02, 1.056, 1.175, 1.002, 0.982, 0.978, 0.947, 1.003, 1.465, 0.985, 1.159, 1.366, 1.309, 0.925, 1.086, 1.212, 1.054, 1.193, 1.079, 1.217, 0.889}, {0.852, 0.735, 1.371, 1.236, 0.926, 1.064, 1.184, 1.323, 0.949, 1.215, 1.224, 1.061, 1.316, 1.186, 1.052, 1.165, 1.316, 0.968, 1.304, 1.076, 1.1, 1.012}, {0.815, 0.716, 1.385, 1.266, 0.974, 1.093, 1.256, 1.337, 0.963, 1.187, 1.235, 1.03, 1.326, 1.161, 1.023, 1.179, 1.292, 0.938, 1.342, 1.09, 1.113, 0.999}, {0.698, 0.662, 1.318, 1.28, 1.093, 1.146, 1.307, 1.367, 0.99, 1.213, 1.217, 0.903, 1.333, 1.086, 0.959, 1.165, 1.204, 0.83, 1.594, 1.174, 1.162, 0.968}, {1.023, 1.002, 0.123, 1.156, 1.03, 1.012, 1.023, 1.031, 1.073, 1.311, 1.08, 1.15, 1.213, 1.156, 1.023, 1.076, 1.142, 1.083, 1.161, 1.076, 1.11, 0.965}, {0.851, 0.768, 1.354, 1.025, 1.181, 1.326, 1.048, 1.226, 0.977, 0.981, 1.237, 1.074, 1.562, 1.694, 1.246, 1.283, 1.105, 0.766, 1.121, 0.878, 1.326, 1.078}, {1.677, 1.722, 0.76, 0.96, 1.2, 1.295, 2.212, 1.296, 0.688, 1.084, 0.872, 0.825, 1.006, 1.216, 1.244, 0.861, 0.977, 0.922, 0.958, 0.926, 1.146, 0.99}};{n, nyear} = Dimensions[R];ER = Mean[Transpose@R]Covariants = Covariance[Transpose[R]];vars = Map[Subscript[x, #]&, {"3m T-bill", "long T-bond", "SP500", "Wiltshire 5000", "Corporate Bond", "NASDQ", "EAFE", "Gold"}];
vars = Map[Subscript[x, #]&, {"3m T-bill", "long T-bond", "SP500", "Wiltshire 5000", "Corporate Bond", "NASDQ", "EAFE", "Gold"}];
FindMinimum[{
vars.Covariants.vars, Total[vars] == 1 && vars.ER ≥ 1.12 && Apply[And, Thread[Greater[vars, 0]]]}, vars]vars = Map[Subscript[x, #]&, {"3m T-bill", "long T-bond", "SP500", "Wiltshire 5000", "Corporate Bond", "NASDQ", "EAFE", "Gold"}];
FindMinimum[{
vars.Covariants.vars, Total[vars] == 1 && vars.ER ≥ 1.10 && Apply[And, Thread[Greater[vars, 0]]]}, vars]内点法の制限
FindMinimumで内点法を実装するには,目的関数と制約条件の第1および第2導関数が必要である.まず記号的導関数が試され,失敗すると有限差分を使って導関数が計算される.関数または制約条件が滑らかでない場合,特に最適点の第1導関数が連続でない場合は,内点法では収束しにくい可能性がある.
FindMinimum[{Abs[x ^ 2 - 3], 0 ≤ x ≤ 5}, {x}, Method -> "InteriorPoint"]
FindMinimum[{Abs[x ^ 2 - 3]}, {x}, Method -> "Newton"]
制約条件付き局所的最適化の数値アルゴリズム
内点法
内点法は障壁関数を使って制約条件と目的関数を組み合せることにより,制約条件付き最適化を解く.特に,一般の制約条件付き最適化問題は,まず次の標準形式
および
に変換される.次に非負の制約条件が,目的関数に障壁項を加えることにより置き換えられる.
例えば,
の勾配が線形非依存であると仮定すると,必要なKKT条件は
であり,
次元
×
である.また
は双対変数である.上記は以下のように書ける:
ここで
は,対角要素が
ならば
の
,あるいは
の対角行列である.また,
は1のベクトルであり
である.変数
を導入すると,上記の系は次のように書き換えられる:
また
は対角要素が
の対角行列である.この非線形系はニュートン法を使って解くことができる.
および
が成り立つとすると,上記の系(1)のヤコビ行列は次の通りである.
で反復的に解くことができる.このとき検索方向
は上記のヤコビ系(2)を解くことで与えられる.
収束を確実なものにするために,成功の測定を行う必要がある.この方法のひとつにメリット関数を使うというものがある.次の拡張ラグランジュ(Lagrangian)メリット関数を使う.
拡張ラグランジュメリット関数
ここで
は障壁パラメータ,
はペナルティパラメータである.行列
が正定値行列ならば,(3)により与えられる検索方向が上記のメリット関数(4)について降順の方向であるか,
がKKT条件(5)を満足するかのいずれかとなることが証明できる.正の制約条件を維持したままで,最初のステップ長ができるだけ1の近くに選ばれて線形探索が検索方向に沿って実行される.次に,Armijo条件がメリット関数
(ここで
)で満足するまで,バックトラック法が使われる.
収束誤差許容値
であり,tolはデフォルトで10-MachinePrecision/3に設定されている.
関連するテクニカルノート
-
▪
- 制約条件付き最適化 ▪
- 線形最適化 ▪
- 大域的非線形数値最適化 ▪
- 厳密な大域的最適化