勾配またはヤコビアンの有効性を確認
最適化ソルバーでは、多くの場合、目的関数と非線形制約関数の 1 次導関数を与えると、計算の速度と信頼性が向上します。ただし、このような導関数のプログラミングではミスが生じやすくなります。導関数の誤りを防ぐには、関数 checkGradients
を使用して、プログラミングした導関数の出力と有限差分近似を照合します。
目的関数の勾配の確認
次の Rastrigin 目的関数について考えます。
次の関数は 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 回行ったところで原則以前と同じ解に辿り着きます。
ベクトル目的関数のヤコビアンの確認
関数 lsqnonlin
、lsqcurvefit
、および fsolve
は、ベクトル値目的関数を使用します。(これらのソルバーが最小化しようとする目的関数は、ベクトル値目的の二乗和です。)"ヤコビアン" は、ベクトル値目的関数の 1 次導関数です。m
個の要素をもつ関数 F(x)
のヤコビアン J
は (ここで、変数 x
には n
個の要素があります)、m
行 n
列の行列です。J(i,j)
は、Fi(x)
の xj
に関する偏導関数です。
たとえば、vecfun
コードは、変数の数が入力データ t
によって異なる 4 つのパラメーター (コード内の a
、b
、c
、および 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
名前と値の引数 IsConstraint
を true
に設定し、点 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.