問題ベースのワークフローへの導関数の供給

導関数を含める理由

この例では、自動導関数が適用されない場合または自動微分の使用では提供されないヘッシアンを含めたい場合に、非線形問題ベースの最適化に導関数情報を含める方法を示します。最適化に勾配またはヘッシアンを含めると、次の利点があります。

より確実な結果が得られます。有限差分のステップでは、目的関数または非線形制約関数が未定義、有限ではない、または複雑である点に到達することがあります。
精度が向上します。解析勾配は、有限差分の推定値より精度が高い場合があります。
収束が速くなります。ヘッシアンを含めることによって収束が早くなる、つまり反復回数が少なくなる可能性があります。
パフォーマンスが改善されます。解析勾配は、特にスパース構造を持つ問題の場合に、有限差分の推定値より早く計算できます。ただし、複雑な式の場合には、解析勾配は計算が遅くなることがあります。

最適化に適用される自動微分

R2020b から、関数 solve がソルバーに勾配を供給するためにサポート対象関数に対して自動微分を使用できるようになりました。これらの自動導関数は、ヘッシアン (2 次導関数) ではなく、勾配 (1 次導関数) にのみ適用されます。

R2022b 以降、自動微分は fcn2optimexpr を使用するかどうかに関係なく、目的関数または制約関数を作成する場合に適用されます。以前は、fcn2optimexpr を使用すると自動微分が適用されませんでした。fcn2optimexpr を使用する必要がある場合に、この例が導関数情報を含める方法を示します。

自動微分を使用せずに問題ベースの最適化で導関数を使用するには、prob2struct を使用して問題を変換してから、結果として得られる目的関数と制約関数を編集します。この例では、自動微分が目的関数の一部として導関数を供給するハイブリッドアプローチを示します。

最適化問題の作成

制御変数 x および y を使って、次の目的関数を使用します。

$\begin{array}{l} fun1 = 100 {(y - x^{2})}^{2} + {(1 - x)}^{2} \\ fun2 = besselj (1, x^{2} + y^{2}) \\ objective = fun1 + fun2 . \end{array}$

x と y の二乗和が 4 以下になるという制約を含めます。

fun2 は最適化式向けにサポートされている関数に基づくものではありません。最適化変数および式でサポートされる演算を参照してください。そのため、fun2 を最適化問題に含めるには、fcn2optimexpr を使用して最適化式に変換しなければなりません。

サポート対象関数に対して AD を使用するには、未サポートの関数 fun2 を使用せずに問題を設定してから、fun2 を含めます。

prob = optimproblem;
x = optimvar('x');
y = optimvar('y');
fun1 = 100*(y - x^2)^2 + (1 - x)^2;
prob.Objective = fun1;
prob.Constraints.cons = x^2 + y^2 <= 4;

ソルバーベース形式への問題の変換

fun2 の導関数を含めるには、まず、prob2struct を使用して、fun2 を使用しない問題を構造体に変換します。

problem = prob2struct(prob,...
    'ObjectiveFunctionName','generatedObjectiveBefore');

変換中に、prob2struct は目的関数と非線形制約関数を表す関数ファイルを作成します。既定では、これらの関数の名前は 'generatedObjective.m' と 'generatedConstraints.m' です。fun2 を使用しない目的関数ファイルは 'generatedObjectiveBefore.m' です。

生成された目的関数と制約関数に勾配が含まれています。

導関数の計算と変数の追跡

fun2 に関連付けられた導関数を計算します。Symbolic Math Toolbox™ ライセンスをお持ちの場合は、関数 gradient (Symbolic Math Toolbox) または jacobian (Symbolic Math Toolbox) を使用して導関数の計算に役立てることができます。詳細については、Symbolic Math Toolbox を使用した勾配とヘッシアンの計算を参照してください。手動で生成した勾配関数を有限差分近似と比較して検証するには、checkGradients を使用します。

ソルバーベースのアプローチには 1 つの制御変数があります。各最適化変数 (この例では x と y) は制御変数の一部です。生成された目的関数ファイル 'generatedObjectiveBefore.m' で、最適化変数と単一の制御変数の間のマッピングを確認できます。ファイルの先頭に次のコード行または同様のコード行が含まれています。

%% Variable indices.
xidx = 1;
yidx = 2;

%% Map solver-based variables to problem-based.
x = inputVariables(xidx);
y = inputVariables(yidx);

このコードで、制御変数の名前は inputVariables です。

または、varindex を使用してマッピングを確認できます。

idx = varindex(prob);
disp(idx.x)

disp(idx.y)

完全な目的関数には、fun2 が含まれています。

fun2 = besselj(1,x^2 + y^2);

標準微積分を使用して、fun2 の勾配である gradfun2 を計算します。

$gradfun2 = [\begin{matrix} 2 x (besselj (0, x^{2} + y^{2}) - besselj (1, x^{2} + y^{2}) / (x^{2} + y^{2})) \\ 2 y (besselj (0, x^{2} + y^{2}) - besselj (1, x^{2} + y^{2}) / (x^{2} + y^{2})) \end{matrix}] .$

目的関数ファイルと制約ファイルの編集

fun2 を含めるように 'generatedObjectiveBefore.m' を編集します。

%% Compute objective function.
arg1 = (y - x.^2);
arg2 = 100;
arg3 = arg1.^2;
arg4 = (1 - x);
obj = ((arg2 .* arg3) + arg4.^2);

ssq = x^2 + y^2;
fun2 = besselj(1,ssq);
obj = obj + fun2;

目的関数ファイルに計算された勾配を含めるには、'generatedObjectiveBefore.m' を編集します。勾配計算を実行しないソフトウェアバージョンを使用している場合は、以下のすべての行を含めてください。ソフトウェアが勾配計算を実行する場合は、fun2 の勾配を計算する太字の行だけを含めてください。

%% Compute objective gradient.
if nargout > 1
    arg5 = 1;
    arg6 = zeros([2, 1]);
    arg6(xidx,:) = (-(arg5.*2.*(arg4(:)))) + ((-((arg5.*arg2(:)).*2.*(arg1(:)))).*2.*(x(:)));
    arg6(yidx,:) = ((arg5.*arg2(:)).*2.*(arg1(:)));
    grad = arg6(:);
    
    arg7 = besselj(0,ssq);
    arg8 = arg7 - fun2/ssq;
    gfun = [2*x*arg8;...
        2*y*arg8];
    
    grad = grad + gfun;
end

非線形制約は x^2 + y^2 <= 4 であることを思い出してください。この制約関数の勾配は 2*[x;y] です。ソフトウェアが制約勾配を計算して、それを生成された制約ファイルに含める場合は、これ以上何もする必要はありません。ソフトウェアが制約勾配を計算しない場合は、以下の行を含めるように 'generatedConstraints.m' を編集することによって、非線形制約の勾配を含めてください。

%% Insert gradient calculation here.
% If you include a gradient, notify the solver by setting the
% SpecifyConstraintGradient option to true.
if nargout > 2
    cineqGrad = 2*[x;y];
    ceqGrad = [];
end

勾配を含む問題と含まない問題の実行

ソルバーベースの方法と問題ベース (勾配なし) の方法の両方を使用して問題を実行し、違いを確認します。導関数情報を使用してソルバーベースの問題を実行するには、問題構造体で適切なオプションを作成します。

options = optimoptions('fmincon','SpecifyObjectiveGradient',true,...
    'SpecifyConstraintGradient',true);
problem.options = options;

非線形問題では、問題構造体の x0 フィールドが空でないことが必要です。

x0 = [1;1];
problem.x0 = x0;

問題構造体に対して fmincon を呼び出します。

[xsolver,fvalsolver,exitflagsolver,outputsolver] = fmincon(problem)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

xsolver =

    1.2494
    1.5617


fvalsolver =

   -0.0038


exitflagsolver =

     1


outputsolver = 

  struct with fields:

         iterations: 15
          funcCount: 32
    constrviolation: 0
           stepsize: 1.5569e-06
          algorithm: 'interior-point'
      firstorderopt: 2.2058e-08
       cgiterations: 7
            message: '↵Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 2.125635e-08,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
       bestfeasible: [1×1 struct]

この解と、導関数情報を使用せずに solve から取得された解を比較します。

init.x = x0(1);
init.y = x0(2);
f2 = @(x,y)besselj(1,x^2 + y^2);
fun2 = fcn2optimexpr(f2,x,y);
prob.Objective = prob.Objective + fun2;
[xproblem,fvalproblem,exitflagproblem,outputproblem] = solve(prob,init,...
    "ConstraintDerivative","finite-differences",...
    "ObjectiveDerivative","finite-differences")

Solving problem using fmincon.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

xproblem = 

  struct with fields:

    x: 1.2494
    y: 1.5617


fvalproblem =

   -0.0038


exitflagproblem = 

    OptimalSolution


outputproblem = 

  struct with fields:

              iterations: 15
               funcCount: 64
         constrviolation: 0
                stepsize: 1.5571e-06
               algorithm: 'interior-point'
           firstorderopt: 6.0139e-07
            cgiterations: 7
                 message: '↵Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 5.795368e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
            bestfeasible: [1×1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

どちらの解でも表示精度に対する関数値は同じ値が報告されており、どちらの解も同じ回数の反復が必要です。ただし、勾配を使用しない解では 64 回の関数評価が必要であるのに対し、勾配を使用する解では関数評価が 32 回しか必要ありません。

ヘッシアンの包含

ヘッシアンを含めるには、すべての関数が最適化式に対してサポートされている場合でも、prob2struct を使用しなければなりません。この例では、fmincon interior-point アルゴリズムにヘッシアンを使用する方法を示します。fminunc trust-region アルゴリズムと fmincon trust-region-reflective アルゴリズムは別々の構文を使用します。fminunc の信頼領域法アルゴリズムまたは fmincon の信頼領域 Reflective 法アルゴリズムのヘッシアンを参照してください。

fmincon の内点法アルゴリズムのヘッシアンで説明されているように、ヘッシアンはラグランジュ関数のヘッシアンです。

\nabla_{x x}^{2} L (x, λ) = \nabla^{2} f (x) + \sum λ_{g, i} \nabla^{2} g_{i} (x) + \sum λ_{h, i} \nabla^{2} h_{i} (x) .

(1)

ヘッシアンを計算するための関数ファイルを作成し、fmincon でそのヘッシアンを使用するための HessianFcn オプションを作成することによって、このヘッシアンを含めます。このようにヘッシアンを作成するには、目的と非線形制約の 2 次導関数を別々に作成します。

この例の目的関数の 2 次導関数行列はやや複雑です。その目的関数リスト hessianfun(x) は、Symbolic Math Toolbox を使用した勾配とヘッシアンの計算で説明されているアプローチと同じアプローチを使用して Symbolic Math Toolbox によって作成されました。

function hf = hessfun(in1)
%HESSFUN
%    HF = HESSFUN(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.6.
%    10-Aug-2020 10:50:44

x = in1(1,:);
y = in1(2,:);
t2 = x.^2;
t4 = y.^2;
t6 = x.*4.0e+2;
t3 = t2.^2;
t5 = t4.^2;
t7 = -t4;
t8 = -t6;
t9 = t2+t4;
t10 = t2.*t4.*2.0;
t11 = besselj(0,t9);
t12 = besselj(1,t9);
t13 = t2+t7;
t14 = 1.0./t9;
t16 = t3+t5+t10-2.0;
t15 = t14.^2;
t17 = t11.*t14.*x.*y.*4.0;
t19 = t11.*t13.*t14.*2.0;
t18 = -t17;
t20 = t12.*t15.*t16.*x.*y.*4.0;
t21 = -t20;
t22 = t8+t18+t21;
hf = reshape([t2.*1.2e+3-t19-y.*4.0e+2-t12.*t15.*...
    (t2.*-3.0+t4+t2.*t5.*2.0+t3.*t4.*4.0+t2.^3.*2.0).*2.0+2.0,...
    t22,t22,...
    t19-t12.*t15.*(t2-t4.*3.0+t2.*t5.*4.0+...
    t3.*t4.*2.0+t4.^3.*2.0).*2.0+2.0e+2],[2,2]);

対照的に、非線形不等式制約のヘッシアンは単純です。単位行列の 2 倍です。

hessianc = 2*eye(2);

ラグランジュ関数のヘッシアンを関数ハンドルとして作成します。

H = @(x,lam)(hessianfun(x) + hessianc*lam.ineqnonlin);

このヘッシアンを使用するためのオプションを作成します。

problem.options.HessianFcn = H;

問題を解いて、反復の回数および関数評価の回数を表示します。解は前回とほぼ同じです。

[xhess,fvalhess,exitflaghess,outputhess] = fmincon(problem);
disp(outputhess.iterations)
disp(outputhess.funcCount)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

    8

    10

今回、fmincon は、15 回実行されていた反復を 8 回しか実行せず、32 回実行されていた関数評価を 10 回しか実行しませんでした。要約すると、解析的ヘッシアン計算を追加することによって、解法プロセスの効率を高めることはできますが、ヘッシアンを計算するための関数の作成が複雑になる可能性があります。

参考

prob2struct | varindex | fcn2optimexpr | checkGradients