ドキュメンテーション

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

Symbolic Math Toolbox による勾配とヘッセ行列の計算

Symbolic Math Toolbox™ のライセンスをおもちの場合は、目的関数と制約関数の解析勾配とヘッセ行列を簡単に計算できます。2 つの関連する Symbolic Math Toolbox 関数があります。

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

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

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

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

z|x||y|(6-60)
x2+y2+(z+1)21.(6-61)

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

 図を生成するコード

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

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

勾配とヘッセ行列を効率よく使用する方法については、勾配とヘッセ行列なしの最適化と比較を参照してください。

変数を作成します。

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

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 = 
  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

線形制約を含める

式 6-60にある線形制約を記述します。

z ≤ –|x| – |y|,

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

xi3xi1xi2 ≤ 0
xi3xi1 + xi2 ≤ 0
xi3 + xi1xi2 ≤ 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 が不等式を表すことがわかります。

A*x
 
ans = 
    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

非線形制約とその勾配とヘッセ行列を作成

式 6-61にある非線形制約

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 のように別個の変数にするよりも、ここで示すように、セル配列でそれらを格納するほうが得策です。

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

目的関数、その勾配、ヘッセ行列の作成

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

energy=i<j1|xixj|.

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

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

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

目的関数と制約用のファイルをすべて生成したら、以下のように適切なラグランジュ乗数で、それらを M ファイル hessfinal.m 内で一緒にします。

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);

最適化の実行

[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','GradObj','on',...
    'GradConstr','on','Hessian','user-supplied',...
    'HessFcn',@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 default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

xfinal =
   -0.0317
    0.0317
   -1.9990
    0.6356
   -0.6356
   -1.4381
    0.0000
   -0.0000
   -0.0000
    0.0000
   -1.0000
   -1.0000
    1.0000
   -0.0000
   -1.0000
   -1.0000
   -0.0000
   -1.0000
    0.6689
    0.6644
   -1.3333
   -0.6667
    0.6667
   -1.3333
    0.0000
    1.0000
   -1.0000
   -0.6644
   -0.6689
   -1.3333

fval =
   34.1365

exitflag =
     1

output = 
         iterations: 19
          funcCount: 28
    constrviolation: 0
           stepsize: 4.0372e-05
          algorithm: 'interior-point'
      firstorderopt: 4.0015e-07
       cgiterations: 55
            message: 'Local minimum found that satisfies the constraints.

Optimizati...'

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

 図を生成するコード

勾配とヘッセ行列なしの最適化と比較

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

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

出力は、fmincon が等しい最小を見つけたものの、それに至るまで、より多くの反復とより多くの関数評価を行ったことを示します。

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.

xfinal2 =
    0.0000
    1.0000
   -1.0000
    0.6689
   -0.6644
   -1.3334
   -0.6644
    0.6689
   -1.3334
    0.0000
   -1.0000
   -1.0000
    0.6357
    0.6357
   -1.4380
   -0.0317
   -0.0317
   -1.9990
    1.0000
    0.0000
   -1.0000
   -1.0000
    0.0000
   -1.0000
    0.0000
    0.0000
   -0.0000
   -0.6667
   -0.6667
   -1.3334

fval2 =
   34.1365

exitflag2 =
     1

output2 = 
         iterations: 77
          funcCount: 2431
    constrviolation: 0
           stepsize: 6.0588e-07
          algorithm: 'interior-point'
      firstorderopt: 2.9894e-06
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.

Optimizati...'

この実行では関数評価 (output2.funcCount) が 2431 回実行されます。これに対し、勾配とヘッセ行列を使用した場合は 28 回 (output.funcCount) になります。

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

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

syms x11 x12 x13 clear 

または、次のコマンドでシンボリック計算エンジンをリセットしなければなりません。

reset(symengine)

シンボリック計算エンジンをリセットした後、clear コマンドか clear variable_list コマンドで MATLAB ワークスペースからすべてのシンボリック変数をクリアしなければなりません。

この情報は役に立ちましたか?