Main Content

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

最適化ソルバーでは、多くの場合、目的関数と非線形制約関数の 1 次導関数を与えると、計算の速度と信頼性が向上します。ただし、このような導関数のプログラミングではミスが生じやすくなります。導関数の誤りを防ぐには、関数 checkGradients を使用して、プログラミングした導関数の出力と有限差分近似を照合します。

目的関数の勾配の確認

次の Rastrigin 目的関数について考えます。

Ras(x)=20+x12+x2210cos(2πx1+2πx2).

次の関数は Rastrigin 関数とその勾配を正しく計算します。

function [f,g] = ras(x)
f = 20 + x(1)^2 + x(2)^2 - 10*cos(2*pi*x(1) + 2*pi*x(2));
if nargout > 1
    g(2) = 2*x(2) + 20*pi*sin(2*pi*x(1) + 2*pi*x(2));
    g(1) = 2*x(1) + 20*pi*sin(2*pi*x(1) + 2*pi*x(2));
end
end

関数 ras でランダムな初期点近傍の勾配が正しく計算されることを確認します。計算の詳細を確認するには、名前と値のパラメーター Display"on" に設定します。

rng default
x0 = randn(1,2);
valid = checkGradients(@ras,x0,Display="on")
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 7.32446e-08.

checkGradients successfully passed.
____________________________________________________________


valid =

  logical

   1

関数 badras は勾配を正しく計算しません。

function [f,g] = badras(x)
f = 20 + x(1)^2 + x(2)^2 - 10*cos(2*pi*x(1) + 2*pi*x(2));
if nargout > 1
    g(2) = 2*x(2) + 40*pi*sin(2*pi*x(1) + 2*pi*x(2));
    g(1) = 2*x(1) + 40*pi*sin(2*pi*x(1) + 2*pi*x(2));
end
end

checkGradients は、関数 badras が勾配を正しく計算しないことを正しく報告します。

valid = checkGradients(@badras,x0,Display="on")
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 0.494224.
  Supplied derivative element (1,1): 92.9841
  Finite-difference derivative element (1,1): 47.0291

checkGradients failed.

Supplied derivative and finite-difference approximation are
not within 'Tolerance' (1e-06).
____________________________________________________________


valid =

  logical

   0

正しい勾配を使用して、ras(x) の局所的最小値を求めます。

rng default
x0 = randn(1,2);
options = optimoptions("fminunc",SpecifyObjectiveGradient=true);
[x,fval,exitflag,output] = fminunc(@ras,x0,options)
Local minimum found.

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

x =

    0.9975    0.9975


fval =

   11.9949


exitflag =

     1


output = 

  struct with fields:

       iterations: 9
        funcCount: 13
         stepsize: 5.7421e-05
     lssteplength: 1
    firstorderopt: 1.2426e-05
        algorithm: 'quasi-newton'
          message: 'Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The first-order optimality measure, 2.482746e-07, is less ↵than options.OptimalityTolerance = 1.000000e-06.'

ソルバーは反復を 9 回行い、関数評価を 13 回だけ行ったところで局所的最小値に到達します。勾配を使用しないと、ソルバーが行う関数評価の回数が増えます。

[x,fval,exitflag,output] = fminunc(@ras,x0) % No options
Local minimum found.
...
output = 

  struct with fields:

       iterations: 9
        funcCount: 39
...

今回、ソルバーは、関数評価を 39 回行ったところで原則以前と同じ解に辿り着きます。

ベクトル目的関数のヤコビアンの確認

関数 lsqnonlinlsqcurvefit、および fsolve は、ベクトル値目的関数を使用します。(これらのソルバーが最小化しようとする目的関数は、ベクトル値目的の二乗和です。)"ヤコビアン" は、ベクトル値目的関数の 1 次導関数です。m 個の要素をもつ関数 F(x) のヤコビアン J は (ここで、変数 x には n 個の要素があります)、mn 列の行列です。J(i,j) は、Fi(x)xj に関する偏導関数です。

たとえば、vecfun コードは、変数の数が入力データ t によって異なる 4 つのパラメーター (コード内の abc、および d) の関数の勾配を計算します。t ベクトルは一連の時間に対応し、各時間は目的関数 F(x,t) の 2 つのエントリになります。

function [F,J] = vecfun(x,t)
t = t(:); % Reshape t to a column vector
a = x(1);
b = x(2);
c = x(3);
d = x(4);
nt = length(t);
F = [a*exp(-b*t) + c,...
    c*exp(-d*t)]; % numel(F) = 2*nt
if nargout > 1
    J = zeros(2*nt,4);
    J(1:nt,1) = exp(-b*t); % First nt components corresponding to a*exp(-b*t) + c
    J(1:nt,2) = (-t.*a.*exp(-b*t));
    J(1:nt,3) = ones(nt,1);
    J(nt+1:2*nt,3) = exp(-d*t); % Second nt components corresponding to c*exp(-d*t)
    J(nt+1:2*nt,4) = (-t.*c.*exp(-d*t));
end
end

vecfun の勾配計算を検証します。

t = [-1/2 1/2 1]; % Three times
x = [1 2 3 4]; % Arbitrary x point
valid = checkGradients(@(x)vecfun(x,t),x,Display="on")
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.51598e-08.

checkGradients successfully passed.
____________________________________________________________


valid =

  logical

   1

一連の時間における関数 vecfun に対応する人為的な応答データを作成します。ヤコビアンを使用して関数のパラメーターをデータに当てはめます。

t = linspace(1,5);
x0 = [1 2 3 4];
rng default
x01 = rand(1,4);
y = vecfun(x0,t);
y = y + 0.1*randn(size(y)); % Add noise to response
options = optimoptions("lsqcurvefit",SpecifyObjectiveGradient=true);
[x,resnorm,~,exitflag,output] = lsqcurvefit(@vecfun,x01,t,y,[],[],options)
Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

x =

    0.9575    1.4570    2.9894    3.7289


resnorm =

    2.2076


exitflag =

     3


output = 

  struct with fields:

      firstorderopt: 1.0824e-04
         iterations: 11
          funcCount: 12
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 0.0111
            message: 'Local minimum possible.↵↵lsqcurvefit stopped because the final change in the sum of squares relative to ↵its initial value is less than the value of the function tolerance.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'
       bestfeasible: []
    constrviolation: []

返された解 [0.9575,1.4570,2.9894,3.7289] は、元のベクトル [1,2,3,4] に近いベクトルです。ソルバーは関数評価を 12 回だけ行います。ヤコビアンを使用しないと、ソルバーが行う評価の回数が増えます。

[x,resnorm,~,exitflag,output] = lsqcurvefit(@vecfun,x01,t,y) % No options
Local minimum possible.
...
output = 

  struct with fields:

      firstorderopt: 1.0825e-04
         iterations: 11
          funcCount: 60
...

今回、ソルバーは、関数評価を 60 回行ったところで原則以前と同じ解に辿り着きます。

有限差分オプションと checkGradients 引数の変更

checkGradients の結果が不正確になる場合があります。

  • 勾配が正しい場合に、checkGradients が誤って無効なチェックを報告する可能性があります。通常これは、関数に比較的大きい 2 次導関数があることで、有限差分推定が不正確になっているために生じます。また、誤った報告は、名前と値の引数 Tolerance の値が小さすぎるために発生することもあります。

  • 勾配が正しくない場合に、checkGradients が誤って有効なチェックを報告する可能性があります。通常この誤った報告は、名前と値の引数 Tolerance の値が大きすぎるため、または結果に偶然一致する導関数が含まれているために発生します。

誤って無効なチェックが報告された場合に勾配をより入念に確認するには、有限差分オプションを変更します。たとえば、checkGradients で、ras 勾配が点 [-2,4] から不正確になっていると誤って報告されたとします。

valid = checkGradients(@ras,[-2,4],Display="on")
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.88497e-06.
  Supplied derivative element (1,2): 6.21515
  Finite-difference derivative element (1,2): 6.21516

checkGradients failed.

Supplied derivative and finite-difference approximation are
not within 'Tolerance' (1e-06).
____________________________________________________________


valid =

  logical

   0

FiniteDifferenceType オプションを "central" に設定し、テストをもう一度実行します。

opts = optimoptions("fmincon",FiniteDifferenceType="central");
valid = checkGradients(@ras,[-2,4],opts,Display="on")
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.10182e-09.

checkGradients successfully passed.
____________________________________________________________


valid =

  logical

   1

通常、中心有限差分法の方がより正確です。この場合、計算された勾配と中心有限差分の推定値の相対差は約 1e-9 です。既定の前進有限差分法での相対差は約 2e-6 です。

許容誤差を既定の 1e-6 から 1e-5 に緩めて、勾配の有効性を確認します。

valid = checkGradients(@ras,[-2,4],Tolerance=1e-5,Display="on")
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.88497e-06.

checkGradients successfully passed.
____________________________________________________________


valid =

  logical

   1

この場合、checkGradients は、勾配が有効であると正しく報告します。許容誤差を緩くしても、関数 badras はチェックにパスしません。

valid = checkGradients(@badras,[-2,4],Tolerance=1e-5,Display="on")
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 0.400823.
  Supplied derivative element (1,2): 4.4368
  Finite-difference derivative element (1,2): 6.21516

checkGradients failed.

Supplied derivative and finite-difference approximation are
not within 'Tolerance' (1e-05).
____________________________________________________________


valid =

  logical

   0

非線形制約の導関数の確認

関数 unitdisk2 は、制約関数とその勾配を正しく計算し、x を半径 1 の円板内に維持します。

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
end

関数 unitdiskb は制約関数の勾配を正しく計算しません。

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

if nargout > 2
    gc = [x(1);x(2)]; % Gradient incorrect: off by a factor of 2
    gceq = [];
end
end

名前と値の引数 IsConstrainttrue に設定し、点 x0 = randn(1,2) における勾配を確認します。

rng default
x0 = randn(1,2);
valid = checkGradients(@unitdisk2,x0,IsConstraint=true,Display="on")
____________________________________________________________

Nonlinear inequality constraint derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.53459e-08.

checkGradients successfully passed.
____________________________________________________________


valid =

  1×2 logical array

   1   1

checkGradients は、unitdisk2 の勾配が有効であると正しく報告します。

valid = checkGradients(@unitdiskb,x0,IsConstraint=true,Display="on")
____________________________________________________________

Nonlinear inequality constraint derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.
  Supplied derivative element (2,1): 1.8324
  Finite-difference derivative element (2,1): 3.66479

checkGradients failed.

Supplied derivative and finite-difference approximation are
not within 'Tolerance' (1e-06).
____________________________________________________________


valid =

  1×2 logical array

   0   1

checkGradients は、unitdiskb の勾配が有効でないと正しく報告します。

スクリプトでの勾配の確認

関数 checkGradients を使用して、有限差分近似が対応する勾配関数と一致しない場合にスクリプトを停止できます。

rng default
x0 = randn(1,2);
opts = optimoptions("fmincon",FiniteDifferenceType="central");
assert(checkGradients(@ras,x0,opts,Display="on"))
assert(all(checkGradients(@unitdisk2,x0,opts,...
    IsConstraint=true,Display="on")))
[x,fval] = fmincon(@ras,x0,[],[],[],[],[],[],@unitdisk2)
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 7.52301e-10.

checkGradients successfully passed.
____________________________________________________________

____________________________________________________________

Nonlinear inequality constraint derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.44672e-11.

checkGradients successfully passed.
____________________________________________________________


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.

x =

    0.4987    0.4987


fval =

   10.4987

正しくない制約関数 unitdiskb を使用して同じスクリプトを実行します。スクリプトが早期に停止することに注意してください。

rng default
x0 = randn(1,2);
opts = optimoptions("fmincon",FiniteDifferenceType="central");
assert(checkGradients(@ras,x0,opts,Display="on"))
assert(all(checkGradients(@unitdiskb,x0,opts,...
    IsConstraint=true,Display="on")))
[x,fval] = fmincon(@ras,x0,[],[],[],[],[],[],@unitdiskb)
____________________________________________________________

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 7.52301e-10.

checkGradients successfully passed.
____________________________________________________________

____________________________________________________________

Nonlinear inequality constraint derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 1.
  Supplied derivative element (2,1): 1.8324
  Finite-difference derivative element (2,1): 3.66479

checkGradients failed.

Supplied derivative and finite-difference approximation are
not within 'Tolerance' (1e-06).
____________________________________________________________

Error using assert
Assertion failed.

参考

関連するトピック