Main Content

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

Symbolic Math Toolbox™ による勾配とヘッシアンの計算

この例では、fmincon と Symbolic Math Toolbox の関数を併用して、非線形最適化問題をより迅速かつロバストに解く方法を説明します。Symbolic Math Toolbox のライセンスをお持ちの場合は、目的関数と制約関数の解析勾配とヘッシアンを簡単に計算できます。関連する Symbolic Math Toolbox の関数は次の 2 つです。

  • jacobian (Symbolic Math Toolbox)は、スカラー関数の勾配を生成し、ベクトル関数の偏導関数の行列を生成します。そのため、たとえば、jacobian を勾配に適用することにより、ヘッセ行列 (目的関数の 2 次導関数) を取得することができます。この例の一部は、目的関数と制約関数のシンボリックな勾配とヘッシアンを生成するためにどのように jacobian を使用するかを示します。

  • matlabFunction (Symbolic Math Toolbox)は、シンボリック表現式の値を計算する無名関数またはファイルのどちらかを生成します。この例は、目的関数と制約関数を評価するファイルとその任意の点における導関数を生成するためにどのように matlabFunction を使用するかを示します。

2 つのツールボックスでは、その関数セットの構文や構造体が異なっています。特に、シンボリック変数は実数または複素数のスカラーですが、Optimization Toolbox™ の関数はベクトル引数を渡します。したがって、目的関数、制約およびそのすべての必須導関数を fminconinterior-point アルゴリズムに適した方式でシンボリックに生成するためには、いくつかの手順を踏む必要があります。

勾配とヘッシアンを効率よく使用する方法については、勾配とヘッシアンなしの最適化と比較を参照してください。導関数情報を使用しないこの問題に対する問題ベースのアプローチについては、静電学における制約付き非線形最適化、問題ベースを参照してください。

問題の説明

10 個の電子を本体に置く静電学問題を考えてみましょう。電子は、本体内に存在する制約に従い、その総位置エネルギーを最小化するように自身を配置します。すべての電子が最小値で本体の境界上にあることは既に知られています。電子には区別がないため、この問題では固有の最小値がありません (1 つの解で電子を並べ替えれば、別の有効な解になります)。この例は、Dolan、Moré、Munson の[58]に基づきます。

この例は、以下の不等式によって定義された本体です。

z-|x|-|y|

x2+y2+(z+1)21

この本体は球の上の角錐のように見えます。

[X,Y] = meshgrid(-1:.01:1);
Z1 = -abs(X) - abs(Y);
Z2 = -1 - sqrt(1 - X.^2 - Y.^2);
Z2 = real(Z2);
W1 = Z1; W2 = Z2;
W1(Z1 < Z2) = nan; % only plot points where Z1 > Z2
W2(Z1 < Z2) = nan; % only plot points where Z1 > Z2
hand = figure; % handle to the figure, since we'll plot more later
set(gcf,'Color','w') % white background
surf(X,Y,W1,'LineStyle','none');
hold on
surf(X,Y,W2,'LineStyle','none');
view(-44,18)

形状の上と下の表面の間にわずかな隙間があります。これは、形状を作成するために使用された一般的なプロット ルーチンの人工物です。このルーチンは、他の面に触れる 1 面上の四角形の領域を消去します。

変数の作成

真にシンボリックな変数 xij、1 ~ 10 の i、および 1 ~ 3 の j から構成される 30 行 1 列のベクトルとして、シンボリックなベクトル x を生成します。これらの変数は電子 i の 3 つの座標を表します。xi1x 座標に、xi2y 座標に、xi3z 座標に対応します。

x = cell(3, 10);
for i = 1:10
    for j = 1:3
        x{j,i} = sprintf('x%d%d',i,j);
    end
end
x = x(:); % now x is a 30-by-1 vector
x = sym(x, 'real');

ベクトル x を検証します。

x
x = 

(x11x12x13x21x22x23x31x32x33x41x42x43x51x52x53x61x62x63x71x72x73x81x82x83x91x92x93x101x102x103)[x11; x12; x13; x21; x22; x23; x31; x32; x33; x41; x42; x43; x51; x52; x53; x61; x62; x63; x71; x72; x73; x81; x82; x83; x91; x92; x93; x101; x102; x103]

線形制約を含める

次の線形制約を記述します。

z-|x|-|y|

これは各電子の 4 つのセットの線形不等式です。

xi3 - xi1 - xi2 ≤ 0
xi3 - xi1 + xi2 ≤ 0
xi3 + xi1 - xi2 ≤ 0
xi3 + xi1 + xi2 ≤ 0

そのため、この問題には全部で 40 の線形不等式があります。

構造体で不等式を記述します。

B = [1,1,1;-1,1,1;1,-1,1;-1,-1,1];

A = zeros(40,30);
for i=1:10
    A(4*i-3:4*i,3*i-2:3*i) = B;
end

b = zeros(40,1);

A*x ≤ b が不等式を表すことがわかります。

disp(A*x)

(x11+x12+x13x12-x11+x13x11-x12+x13x13-x12-x11x21+x22+x23x22-x21+x23x21-x22+x23x23-x22-x21x31+x32+x33x32-x31+x33x31-x32+x33x33-x32-x31x41+x42+x43x42-x41+x43x41-x42+x43x43-x42-x41x51+x52+x53x52-x51+x53x51-x52+x53x53-x52-x51x61+x62+x63x62-x61+x63x61-x62+x63x63-x62-x61x71+x72+x73x72-x71+x73x71-x72+x73x73-x72-x71x81+x82+x83x82-x81+x83x81-x82+x83x83-x82-x81x91+x92+x93x92-x91+x93x91-x92+x93x93-x92-x91x101+x102+x103x102-x101+x103x101-x102+x103x103-x102-x101)[x11 + x12 + x13; x12 - x11 + x13; x11 - x12 + x13; x13 - x12 - x11; x21 + x22 + x23; x22 - x21 + x23; x21 - x22 + x23; x23 - x22 - x21; x31 + x32 + x33; x32 - x31 + x33; x31 - x32 + x33; x33 - x32 - x31; x41 + x42 + x43; x42 - x41 + x43; x41 - x42 + x43; x43 - x42 - x41; x51 + x52 + x53; x52 - x51 + x53; x51 - x52 + x53; x53 - x52 - x51; x61 + x62 + x63; x62 - x61 + x63; x61 - x62 + x63; x63 - x62 - x61; x71 + x72 + x73; x72 - x71 + x73; x71 - x72 + x73; x73 - x72 - x71; x81 + x82 + x83; x82 - x81 + x83; x81 - x82 + x83; x83 - x82 - x81; x91 + x92 + x93; x92 - x91 + x93; x91 - x92 + x93; x93 - x92 - x91; x101 + x102 + x103; x102 - x101 + x103; x101 - x102 + x103; x103 - x102 - x101]

非線形制約とその勾配とヘッシアンを作成

次の非線形制約

x2+y2+(z+1)21

は構造体でもあります。制約とその勾配とヘッシアンを以下のように作成します。

c = sym(zeros(1,10));
i = 1:10;
c = (x(3*i-2).^2 + x(3*i-1).^2 + (x(3*i)+1).^2 - 1).';

gradc = jacobian(c,x).'; % .' performs transpose

hessc = cell(1, 10);
for i = 1:10
    hessc{i} = jacobian(gradc(:,i),x);
end

制約ベクトル c は行ベクトルです。c(i) の勾配は、行列 gradci 列目に示されます。非線形制約で説明するように、これは正しい形式です。

ヘッセ行列 hessc{1}, ..., hessc{10} は、平方で対称的です。hessc1, ..., hesssc10 のように別個の変数にするよりも、ここで示すように、cell 配列でそれらを格納するほうが得策です。

.' 構文を使用して、転置します。構文 ' は、異なるシンボリック導関数をもつ共役転置を意味します。

目的関数、その勾配、ヘッシアンの作成

目的関数 (位置エネルギー) は、各電子ペアの間の距離の逆数の合計です。

energy=i<j1|xi-xj|

距離はベクトルの成分にある差分の平方の合計の平方根です。

エネルギー、その勾配とヘッシアンを以下のように計算します。

energy = sym(0);
for i = 1:3:25
    for j = i+3:3:28
        dist = x(i:i+2) - x(j:j+2);
        energy = energy + 1/sqrt(dist.'*dist);
    end
end

gradenergy = jacobian(energy,x).';

hessenergy = jacobian(gradenergy,x);

目的関数ファイルの作成

目的関数は 2 つの出力である energygradenergy をもたなければなりません。matlabFunction を呼び出す際に、両方の関数を 1 つのベクトルに入れます。matlabFunction が生成する部分表現の数が減り、呼び出す関数 (この例では fmincon) が両方の出力を要求する時にのみ、勾配が返されます。この例は、結果のファイルを現在のフォルダーに置いていることを示します。もちろん、フォルダーが MATLAB パス上にあれば、任意の位置にそれらを置くことができます。

currdir = [pwd filesep]; % You may need to use currdir = pwd 
filename = [currdir,'demoenergy.m'];
matlabFunction(energy,gradenergy,'file',filename,'vars',{x});

この構文は、matlabFunction に、最初の出力として energy を、2 番目の出力として gradenergy を返させます。また、入力 x11, ..., x103 のリストの代わりに、単一の入力ベクトル {x} を取ります。

結果ファイル demoenergy.m は、以下の行かそれに似た行を部分的に含みます。

function [energy,gradenergy] = demoenergy(in1)
%DEMOENERGY
%   [ENERGY,GRADENERGY] = DEMOENERGY(IN1)
...
x101 = in1(28,:);
...
energy = 1./t140.^(1./2) + ...;
if nargout > 1
    ...
    gradenergy = [(t174.*(t185 - 2.*x11))./2 - ...];
end

この関数は勾配をもつ目的関数のための正しい形式をもちます (スカラー目的関数の記述を参照)。

制約関数ファイルの作成

非線形制約関数を生成し、それを正しい形式にします。

filename = [currdir,'democonstr.m'];
matlabFunction(c,[],gradc,[],'file',filename,'vars',{x},...
    'outputs',{'c','ceq','gradc','gradceq'});

結果ファイル democonstr.m は、以下の行かそれに似た行を部分的に含みます。

function [c,ceq,gradc,gradceq] = democonstr(in1)
%DEMOCONSTR
%    [C,CEQ,GRADC,GRADCEQ] = DEMOCONSTR(IN1)
...
x101 = in1(28,:);
...
c = [t417.^2 + ...];
if nargout > 1
    ceq = [];
end
if nargout > 2
    gradc = [2.*x11,...];
end
if nargout > 3
    gradceq = [];
end

この関数は勾配をもつ制約関数のための正しい形式をもちます (非線形制約を参照)。

ヘッセ ファイルの生成

問題用のラグランジュ関数のヘッシアンを生成するには、まず、エネルギー ヘッシアンと制約ヘッシアンのファイルを生成します。

目的関数 hessenergy のヘッシアンは、非常に大きいシンボリック表現式です。size(char(hessenergy)) の評価が示すように、150,000 を超えるシンボルが含まれています。そのため、matlabFunction(hessenergy) を実行するために長い時間がかかります。

ファイルの hessenergy.m を生成するには、以下の 2 つの行を実行します。

filename = [currdir,'hessenergy.m'];
matlabFunction(hessenergy,'file',filename,'vars',{x});

その反対に制約関数のヘッシアンは小さく、迅速に計算します。

for i = 1:10
    ii = num2str(i);
    thename = ['hessc',ii];
    filename = [currdir,thename,'.m'];
    matlabFunction(hessc{i},'file',filename,'vars',{x});
end

目的関数と制約用のファイルをすべて生成したら、以下のように適切なラグランジュ乗数で、それらをファイル hessfinal.m (この例の終わりにコードを掲載) 内で一緒にします。

最適化の実行

[0,0,–1] でセンタリングされた半径 1/2 の球上で無作為に分布された電子で最適化を開始します。

rng default % for reproducibility
Xinitial = randn(3,10); % columns are normal 3-D vectors
for j=1:10
    Xinitial(:,j) = Xinitial(:,j)/norm(Xinitial(:,j))/2;
    % this normalizes to a 1/2-sphere
end
Xinitial(3,:) = Xinitial(3,:) - 1; % center at [0,0,-1]
Xinitial = Xinitial(:); % Convert to a column vector

オプションが内点法アルゴリズムで勾配とヘッシアンを使用するように設定します。

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

fmincon を呼び出します。

[xfinal fval exitflag output] = fmincon(@demoenergy,Xinitial,...
    A,b,[],[],[],[],@democonstr,options);
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>

目的関数値、終了フラグ、反復の回数、関数評価の回数を表示します。

disp(fval)
   34.1365
disp(exitflag)
     1
disp(output.iterations)
    19
disp(output.funcCount)
    28

電子の初期の位置が無作為でも、最終的な位置はほぼ対称的です。

for i = 1:10
    plot3(xfinal(3*i-2),xfinal(3*i-1),xfinal(3*i),'r.','MarkerSize',25);
end

rotate3d
figure(hand)

勾配とヘッシアンなしの最適化と比較

勾配とヘッシアンを使用すると、最適化をより速くより正確に実行できます。勾配やヘッセ情報を使用しない最適化と比較するために、勾配とヘッシアンを使用しないようにオプションを設定します。

options = optimoptions(@fmincon,'Algorithm','interior-point',...
    'Display','final');
[xfinal2 fval2 exitflag2 output2] = fmincon(@demoenergy,Xinitial,...
    A,b,[],[],[],[],@democonstr,options);
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>

前回と同じデータを検証します。

disp(fval2)
   34.1365
disp(exitflag2)
     1
disp(output2.iterations)
    78
disp(output2.funcCount)
        2463

導関数情報を使用する場合と使用しない場合で、反復の回数および関数評価の回数を比較します。

tbl = table([output.iterations;output2.iterations],[output.funcCount;output2.funcCount],...
    'RowNames',{'With Derivatives','Without Derivatives'},'VariableNames',{'Iterations','Fevals'})
tbl=2×2 table
                           Iterations    Fevals
                           __________    ______

    With Derivatives           19           28 
    Without Derivatives        78         2463 

シンボリックな変数前提条件をクリア

この例のシンボリックな変数は、シンボリックなエンジン ワークスペース内でそれらが実数であるという前提条件をもっています。この想定をシンボリック エンジン ワークスペースから取り除くには、これらの変数を削除するだけでは不十分です。syms を使用して、変数の仮定をクリアします。

syms x

仮定が空になっていることを確認します。

assumptions(x)
 
ans =
 
Empty sym: 1-by-0
 

補助関数

次のコードは関数 hessfinal を作成します。

function H = hessfinal(X,lambda)
%
% Call the function hessenergy to start
H = hessenergy(X);

% Add the Lagrange multipliers * the constraint Hessians
H = H + hessc1(X) * lambda.ineqnonlin(1);
H = H + hessc2(X) * lambda.ineqnonlin(2);
H = H + hessc3(X) * lambda.ineqnonlin(3);
H = H + hessc4(X) * lambda.ineqnonlin(4);
H = H + hessc5(X) * lambda.ineqnonlin(5);
H = H + hessc6(X) * lambda.ineqnonlin(6);
H = H + hessc7(X) * lambda.ineqnonlin(7);
H = H + hessc8(X) * lambda.ineqnonlin(8);
H = H + hessc9(X) * lambda.ineqnonlin(9);
H = H + hessc10(X) * lambda.ineqnonlin(10);
end

関連するトピック