スカラー目的関数の記述

関数ファイル

スカラー目的関数ファイルは、x などの 1 つの入力を受け入れ、f などの 1 つの実数スカラー出力を返します。入力 x には、スカラー、ベクトルまたは行列を指定できます。関数ファイルでは、より多くの出力を返すことができます (勾配とヘッシアンを含めるを参照)。

たとえば、目的が次のような 3 つの変数 x、y、および z の関数だとします。

f(x) = 3*(x – y)⁴ + 4*(x + z)² / (1 + x² + y² + z²) + cosh(x – 1) + tanh(y + z)

1. この関数を、ベクトル xin = [x;y;z] を受け入れ、f を返すファイルとして記述します。

   function f = myObjective(xin)
   f = 3*(xin(1)-xin(2))^4 + 4*(xin(1)+xin(3))^2/(1+norm(xin)^2) ...
       + cosh(xin(1)-1) + tanh(xin(2)+xin(3));

2. これを myObjective.m という名前のファイルとして、MATLAB^® パス上のフォルダーに保存します。

3. 関数が正しく評価されることをチェックします。

   myObjective([1;2;3])
   
   ans =
       9.2666

追加パラメーターを含める方法は 追加パラメーターの受け渡し を参照してください。さらに複雑な関数ファイルの例は、勾配およびヘッセ スパース パターンを使った最小化または範囲制約および範囲付きの前提条件子を使った最小化を参照してください。

ローカル関数と入れ子関数

関数は、他のファイルの内部にローカル関数または入れ子関数として存在できます。ローカル関数や入れ子関数を使用すると、保存するファイルの数を減らせます。また、入れ子関数を使用すると、入れ子関数に示すように追加のパラメーターにアクセスできるようになります。

たとえば、関数ファイルに記述されている目的関数 myObjective.m を、非線形制約に記述されている ellipseparabola.m の制約に基づいて最小化するとします。この場合、2 つのファイル myObjective.m と ellipseparabola.m を作成する代わりに、両方の関数をローカル関数として含む 1 つのファイルを作成できます。

function [x fval] = callObjConstr(x0,options)
% Using a local function for just one file

if nargin < 2
    options = optimoptions('fmincon','Algorithm','interior-point');
end

[x fval] = fmincon(@myObjective,x0,[],[],[],[],[],[], ...
    @ellipseparabola,options);

function f = myObjective(xin)
f = 3*(xin(1)-xin(2))^4 + 4*(xin(1)+xin(3))^2/(1+sum(xin.^2)) ...
    + cosh(xin(1)-1) + tanh(xin(2)+xin(3));

function [c,ceq] = ellipseparabola(x)
c(1) = (x(1)^2)/9 + (x(2)^2)/4 - 1;
c(2) = x(1)^2 - x(2) - 1;
ceq = [];

点 [1;1;1] から開始して、制約付きの最小化を解きます。

[x fval] = callObjConstr(ones(3,1))

Local minimum found that satisfies the constraints.

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

x =
    1.1835
    0.8345
   -1.6439

fval =
    0.5383

無名関数の目的関数

無名関数は、単純な目的関数の記述に使用します。無名関数の詳細については、無名関数とはを参照してください。Rosenbrock 関数はシンプルな関数であり、無名関数として記述するのに適しています。

anonrosen = @(x)(100*(x(2) - x(1)^2)^2 + (1-x(1))^2);

anonrosen が [-1 2] の点で正しく評価されることをチェックします。

anonrosen([-1 2])

ans =
   104

anonrosen を fminunc で最小化すると以下の結果になります。

options = optimoptions(@fminunc,'Algorithm','quasi-newton');
[x fval] = fminunc(anonrosen,[-1;2],options)

Local minimum found.

Optimization completed because the size of the gradient
is less than the default value of the function tolerance.

x =
    1.0000
    1.0000

fval =
  1.2266e-10

勾配とヘッシアンを含める

• ソルバーへの導関数の指定

• 勾配を含める方法

• ヘッシアンを含める

• 導関数を含める利点

• 内点法 fmincon に対する入力ヘッセ近似の選択

ソルバーへの導関数の指定

fmincon と fminunc に対しては、目的関数に勾配を含めることができます。一般的に、ソルバーに勾配を含めるとロバスト性が向上し、場合によっては速度もやや速くなります。詳細は、導関数を含める利点を参照してください。2 次導関数 (ヘッシアン) も含めるには、ヘッシアンを含めるを参照してください。

以下の表に、どのアルゴリズムで勾配とヘッシアンを使用できるかを示します。

ソルバー	アルゴリズム	勾配	ヘッシアン
fmincon	active-set	オプション	なし
	interior-point	オプション	オプション (fmincon の内点法アルゴリズムのヘッシアンを参照)
	sqp	オプション	なし
	trust-region-reflective	必須	オプション (fminunc の信頼領域法アルゴリズムまたは fmincon の信頼領域 Reflective 法アルゴリズムのヘッシアンを参照)
fminunc	quasi-newton	オプション	なし
	trust-region	必須	オプション (fminunc の信頼領域法アルゴリズムまたは fmincon の信頼領域 Reflective 法アルゴリズムのヘッシアンを参照)

勾配を含める方法

1. 以下の出力を返すコードを記述します。

   • 最初の出力として目的関数 (スカラー)

   • 2 番目の出力として勾配 (ベクトル)

2. optimoptions を使用して SpecifyObjectiveGradient オプションを true に設定します。適切な場合は、SpecifyConstraintGradient オプションも true に設定します。

3. オプションで、勾配関数が有限差分近似と一致するかどうかを確認します。詳細は、勾配またはヤコビアンの有効性を確認を参照してください。

ヒント

最大の柔軟性を得るには、"条件付き" コードを記述します。条件付きとは、以下の例に示すように、関数出力の数が異なることを意味します。条件付きコードでは、SpecifyObjectiveGradient オプションの値によってエラーが発生することはありません。条件なしのコードでは、オプションを適切に設定しなければなりません。

たとえば、次の Rosenbrock 関数について考察します。

$$f(x)=100{\left(x_2-x_1^2\right)}^2+{(1-x_1)}^2,$$

この関数は、[最適化] ライブ エディター タスクまたはソルバーを使用した制約付き非線形問題で詳細が述べられています。f(x) の勾配は、次のようになります。

$$\nabla f(x)=\left[\begin{array}{c}-400\left(x_2-x_1^2\right)x_1-2\left(1-x_1\right)\\ 200\left(x_2-x_1^2\right)\end{array}\right],$$

rosentwo は条件付き関数で、ソルバーで必要なものを返します。

function [f,g] = rosentwo(x)
% Calculate objective f
f = 100*(x(2) - x(1)^2)^2 + (1-x(1))^2;

if nargout > 1 % gradient required
    g = [-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1));
        200*(x(2)-x(1)^2)];
    
end

nargout は呼び出し関数が指定する引数の数をチェックします。詳細は、関数の引数の数の確認を参照してください。

制約なしの最適化のために設計された fminunc ソルバーは Rosenbrock 関数を最小化します。options を設定して、fminunc が勾配とヘッシアンを使用するようにします。

options = optimoptions(@fminunc,'Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true);

fminunc を [-1;2] から開始して起動します。

[x fval] = fminunc(@rosentwo,[-1;2],options)
Local minimum found.

Optimization completed because the size of the gradient
is less than the default value of the function tolerance.

x =
    1.0000
    1.0000

fval =
  1.9886e-17

Symbolic Math Toolbox™ のライセンスをおもちの場合は、勾配とヘッシアンを自動的に計算することができます (Symbolic Math Toolbox を使用した勾配とヘッシアンの計算を参照)。

ヘッシアンを含める

fmincon の 'trust-region-reflective' アルゴリズムと 'interior-point' アルゴリズム、および fminunc の 'trust-region' アルゴリズムには、2 次導関数を含めることができます。ヘッセ情報を含める方法は、情報のタイプとアルゴリズムに応じていくつかあります。

ヘッシアンを含めるには、勾配も含めなければなりません (SpecifyObjectiveGradient を true に設定し、可能な場合は SpecifyConstraintGradient を true 設定します)。

• fminunc の信頼領域法アルゴリズムまたは fmincon の信頼領域 Reflective 法アルゴリズムのヘッシアン

• fmincon の内点法アルゴリズムのヘッシアン

• ヘッセ乗算関数

fminunc の信頼領域法アルゴリズムまたは fmincon の信頼領域 Reflective 法アルゴリズムのヘッシアン.  これらのアルゴリズムには、制約がないか、範囲制約と線形等式制約のいずれかのみがあります。そのため、ヘッシアンは目的関数の 2 次導関数の行列です。

ヘッセ行列を目的関数の 3 番目の出力として含めます。たとえば、Rosenbrock 関数のヘッシアン H(x) は次のとおりです (勾配を含める方法を参照)。

$$H(x)=\left[\begin{array}{cc}1200x_1^2-400x_2+2& -400x_1\\ -400x_1& 200\end{array}\right].$$

このヘッシアンを目的関数に含めます。

function [f, g, H] = rosenboth(x)
% Calculate objective f
f = 100*(x(2) - x(1)^2)^2 + (1-x(1))^2;

if nargout > 1 % gradient required
    g = [-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1));
        200*(x(2)-x(1)^2)];
    
    if nargout > 2 % Hessian required
        H = [1200*x(1)^2-400*x(2)+2, -400*x(1);
            -400*x(1), 200];  
    end

end

HessianFcn を 'objective' に設定します。たとえば、

options = optimoptions('fminunc','Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true,'HessianFcn','objective');

fmincon の内点法アルゴリズムのヘッシアン.  このヘッシアンはラグランジュ関数のヘッシアンです。ここで、ラグランジュ関数 L(x,λ) は次のとおりです。

$$L(x,\lambda )=f(x)+{\displaystyle \sum {\lambda }_{g,i}{g}_i(x)}+{\displaystyle \sum {\lambda }_{h,i}{h}_i(x)}.$$

g と h はそれぞれ不等式制約と等式制約 (つまり、範囲、線形、非線形の制約) を表すベクトル関数であり、最小化問題は以下のようになります。

$$\underset{x}{\min }f(x)\text{ subject to }g(x)\le 0,h(x)=0.$$

詳細は、制約付きの最適性の理論を参照してください。ラグランジュ関数のヘッシアンは以下のようになります。

$${\nabla }_{xx}^{2}L(x,\lambda )={\nabla }^{2}f(x)+{\displaystyle \sum {\lambda }_{g,i}{\nabla }^2{g}_i(x)}+{\displaystyle \sum {\lambda }_{h,i}{\nabla }^2{h}_i(x)}.$$

(1)

ヘッシアンを含めるには、次の構文で関数を記述します。

hessian = hessianfcn(x,lambda)

hessianは、n 行 n 列の行列でなければなりません。ここで、n は変数の数です。hessian が大きく、比較的非零の項目が少ない場合は、hessian をスパース行列として使用することにより、実行時間とメモリを節約します。lambda は、非線形制約に関連するラグランジュ乗数ベクトルをもつ構造体です。

lambda.ineqnonlin
lambda.eqnonlin

fmincon は構造体 lambda を計算し、その結果をヘッセ関数に渡します。hessianfcn は式 1 の合計を計算しなければなりません。次のオプションを設定して、ヘッシアンを渡します。

options = optimoptions('fmincon','Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
    'HessianFcn',@hessianfcn);

たとえば、Rosenbrock 関数のヘッシアンを単位円板 $$x_1^2+x_2^2\le 1$$ の制約付きで含めるには、制約関数 $$g(x)=x_1^2+x_2^2-1\le 0$$ が勾配と 2 次導関数行列をもつことに注目します。

$$\begin{array}{c}\nabla g(x)=\left[\begin{array}{c}2x_1\\ 2x_2\end{array}\right]\\ H_g(x)=\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right].\end{array}$$

ヘッセ関数を以下のように記述します。

function Hout = hessianfcn(x,lambda)
% Hessian of objective
H = [1200*x(1)^2-400*x(2)+2, -400*x(1);
            -400*x(1), 200];
% Hessian of nonlinear inequality constraint
Hg = 2*eye(2);
Hout = H + lambda.ineqnonlin*Hg;

MATLAB パス上に hessianfcn を保存します。この例の場合、勾配を含む制約関数は以下のようになります。

function [c,ceq,gc,gceq] = unitdisk2(x)
c = x(1)^2 + x(2)^2 - 1;
ceq = [ ];

if nargout > 2
    gc = [2*x(1);2*x(2)];
    gceq = [];
end

勾配とヘッシアンを含む問題を解きます。

fun = @rosenboth;
nonlcon = @unitdisk2;
x0 = [-1;2];
options = optimoptions('fmincon','Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
    'HessianFcn',@hessianfcn);
[x,fval,exitflag,output] = fmincon(fun,x0,[],[],[],[],[],[],@unitdisk2,options);

内点法のヘッシアンを使用するその他の例は、解析的ヘッシアンを使用した fmincon の内点法アルゴリズムとSymbolic Math Toolbox を使用した勾配とヘッシアンの計算を参照してください。

ヘッセ乗算関数.  fmincon の interior-point アルゴリズムと trust-region-reflective アルゴリズムの両方で、完全なヘッセ関数の代わりに、ヘッセ乗算関数を指定できます。この関数はヘッシアンを直接計算せずに、ヘッシアンとベクトルとの乗算結果を与えます。これによってメモリを節約できます。ヘッセ乗算関数を機能させるには、SubproblemAlgorithm オプションが 'cg' でなければなりません。trust-region-reflective では、これが既定の設定です。

2 つのアルゴリズムの構文は次のように異なります。

• interior-point アルゴリズムでは、構文は次のようになります。

  W = HessMultFcn(x,lambda,v);

  結果 W は積 H*v になります。ここで H は x でのラグランジュ関数のヘッシアンであり (式 1を参照)、lambda は fmincon によって計算されるラグランジュ乗数です。また、v はサイズが n 行 1 列のベクトルです。以下のようにオプションを設定します。

  options = optimoptions('fmincon','Algorithm','interior-point','SpecifyObjectiveGradient',true,... 
      'SpecifyConstraintGradient',true,'SubproblemAlgorithm','cg','HessianMultiplyFcn',@HessMultFcn);

  n 行 1 列のベクトルを返す関数 HessMultFcn を指定します。ここで、n は x の次元数です。HessianMultiplyFcn オプションはヘッシアンを計算せずにヘッシアンとベクトルの乗算結果を渡すことができます。

• trust-region-reflective アルゴリズムには lambda が含まれません。

  W = HessMultFcn(H,v);

  結果 W は H*v に等しくなります (W = H*v)。fmincon は目的関数の 3 番目の出力で返される値として H を渡します (fminunc の信頼領域法アルゴリズムまたは fmincon の信頼領域 Reflective 法アルゴリズムのヘッシアンを参照)。また fmincon は、n 行をもつベクトルまたは 列 v を渡します。v の列数は可変であるため、任意の列数を受け入れるよう HessMultFcn を記述します。H がヘッシアンである必要はなく、W = H*v を計算できる任意のものを使用できます。

  以下のようにオプションを設定します。

  options = optimoptions('fmincon','Algorithm','trust-region-reflective',... 
      'SpecifyObjectiveGradient',true,'HessianMultiplyFcn',@HessMultFcn);

  ヘッセ乗算関数を trust-region-reflective アルゴリズムとともに使用する例は、密に構造化されたヘッシアンと線形等式を使用した最小化を参照してください。

導関数を含める利点

勾配を指定しない場合、ソルバーは有限差分法により勾配を推定します。勾配を指定すると、ソルバーでこの有限差分推定を行う必要がないため、時間を節約でき、より正確になります。ただし、導関数が複雑な場合は、有限差分推定の方が高速になることがあります。また、ソルバーでは近似ヘッシアンが使用され、真のヘッシアンとは大きく異なる場合があります。ヘッシアンを指定すると、少ない反復数で解が得られます。たとえば、Symbolic Math Toolbox を使用した勾配とヘッシアンの計算の結論を参照してください。

制約問題の場合は、勾配を与えることで別の利点があります。ソルバーは、x それ自体は実行可能な x に対し、x の周囲の有限差分は常に実行不可能な点へ導くような点 x に到達することがあります。また、目的関数は実行不可能な点において、複素数出力、Inf、NaN、またはエラーを返すとします。この場合、ソルバーは失敗するか途中で停止する可能性があります。勾配を与えることで、ソルバーは続行できます。この利点を得るには、さらに非線形制約関数の勾配を含め、SpecifyConstraintGradient オプションを true に設定しなければならない場合があります。詳細は、非線形制約を参照してください。

内点法 fmincon に対する入力ヘッセ近似の選択

fmincon の interior-point アルゴリズムには、入力ヘッセ近似を選択する場合に多くのオプションがあります。構文の詳細は、入力としてのヘッシアン を参照してください。次の表にオプションと各オプション間の相対的な特性をまとめます。

ヘッシアン	相対的なメモリ使用量	相対的な効率性
'bfgs' (既定の設定)	高 (大規模問題の場合)	高
'lbfgs'	低から中	中
'fin-diff-grads'	低	中
'HessianMultiplyFcn'	低 (コードにより変化の可能性)	中
'HessianFcn'	可変 (コードによる)	高 (コードによる)

以下の場合を除き、既定の 'bfgs' ヘッシアンを使用します。

• メモリ不足 — 'bfgs' の代わりに 'lbfgs' を試します。独自の勾配を提供できる場合は 'fin-diff-grads' を使用し、SpecifyObjectiveGradient オプションと SpecifyConstraintGradient オプションを true に設定します。

• 効率性を向上させたい場合 — 独自の勾配とヘッシアンを提供します。ヘッシアンを含める、解析的ヘッシアンを使用した fmincon の内点法アルゴリズム、およびSymbolic Math Toolbox を使用した勾配とヘッシアンの計算を参照してください。

'lbfgs' に中程度の効率性しかないのは、2 つの理由があります。1 つの原因は、Sherman-Morrison の更新に比較的多くのリソースが使用さるためです。また、その結果として生じる反復ステップは 'lbfgs' のメモリが限られているためにやや不正確になります。

'fin-diff-grads' と HessianMultiplyFcn に中程度の効率性しかないのは、共役勾配のアプローチを使用するためです。目的関数のヘッシアンは正確に推定されますが、最も正確な反復ステップが生成されるわけではありません。詳細は、fmincon の内点法アルゴリズムと、式 38 を解決する LDL アプローチと共役勾配アプローチについての説明を参照してください。

関連するトピック

• 勾配またはヤコビアンの有効性を確認