勾配およびヘッシアンを使った最小化
この例では、明示的な三重対角ヘッセ行列 を使用して非線形最小化問題を解く方法を説明します。この問題は以下を最小化する を見つけるものです。
ここで = 1000 です。
補助関数 brownfgh
(この例の終わりに掲載) が、、その勾配 、およびそのヘッシアン を計算します。fminunc
ソルバーが導関数情報を使用するよう指定するには、optimoptions
を使用して SpecifyObjectiveGradient
オプションと HessianFcn
オプションを設定します。fminunc
でヘッシアンを使用するには、'trust-region'
アルゴリズムを使用しなければなりません。
options = optimoptions(@fminunc,'Algorithm','trust-region',... 'SpecifyObjectiveGradient',true,'HessianFcn','objective');
パラメーター n
を 1000 に設定し、初期点 xstart
として奇数の要素に対して -1 を設定し、偶数の要素に対して +1 を設定します。
n = 1000; xstart = -ones(n,1); xstart(2:2:n) = 1;
の最小値を求めます。
[x,fval,exitflag,output] = fminunc(@brownfgh,xstart,options);
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
解法および解法プロセスを検証します。
disp(fval)
2.8709e-17
disp(exitflag)
1
disp(output)
iterations: 7 funcCount: 8 stepsize: 0.0039 cgiterations: 7 firstorderopt: 4.7948e-10 algorithm: 'trust-region' message: 'Local minimum found....' constrviolation: []
関数 は、二乗の和であるため、非負となります。解 fval
はほぼゼロであるため、最小値であることは明らかです。終了フラグが 1
であることも fminunc
が解を見つけたことを示しています。output
構造体は、fminunc
がわずか 7 回の反復で解に到達したことを示しています。
解の最大要素と最小要素を表示します。
disp(max(x))
1.1987e-10
disp(min(x))
-1.1987e-10
この解は、すべての要素が x = 0
となる点に極めて近いものとなります。
補助関数
次のコードは、補助関数 brownfgh
を作成します。
function [f,g,H] = brownfgh(x) %BROWNFGH Nonlinear minimization problem (function, its gradients % and Hessian) % Documentation example % Copyright 1990-2008 The MathWorks, Inc. % Evaluate the function. n=length(x); y=zeros(n,1); i=1:(n-1); y(i)=(x(i).^2).^(x(i+1).^2+1)+(x(i+1).^2).^(x(i).^2+1); f=sum(y); % % Evaluate the gradient. if nargout > 1 i=1:(n-1); g = zeros(n,1); g(i)= 2*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2))+... 2*x(i).*((x(i+1).^2).^(x(i).^2+1)).*log(x(i+1).^2); g(i+1)=g(i+1)+... 2*x(i+1).*((x(i).^2).^(x(i+1).^2+1)).*log(x(i).^2)+... 2*(x(i).^2+1).*x(i+1).*((x(i+1).^2).^(x(i).^2)); end % % Evaluate the (sparse, symmetric) Hessian matrix if nargout > 2 v=zeros(n,1); i=1:(n-1); v(i)=2*(x(i+1).^2+1).*((x(i).^2).^(x(i+1).^2))+... 4*(x(i+1).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i).^2).^((x(i+1).^2)-1))+... 2*((x(i+1).^2).^(x(i).^2+1)).*(log(x(i+1).^2)); v(i)=v(i)+4*(x(i).^2).*((x(i+1).^2).^(x(i).^2+1)).*((log(x(i+1).^2)).^2); v(i+1)=v(i+1)+... 2*(x(i).^2).^(x(i+1).^2+1).*(log(x(i).^2))+... 4*(x(i+1).^2).*((x(i).^2).^(x(i+1).^2+1)).*((log(x(i).^2)).^2)+... 2*(x(i).^2+1).*((x(i+1).^2).^(x(i).^2)); v(i+1)=v(i+1)+4*(x(i).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i+1).^2).^(x(i).^2-1)); v0=v; v=zeros(n-1,1); v(i)=4*x(i+1).*x(i).*((x(i).^2).^(x(i+1).^2))+... 4*x(i+1).*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2)).*log(x(i).^2); v(i)=v(i)+ 4*x(i+1).*x(i).*((x(i+1).^2).^(x(i).^2)).*log(x(i+1).^2); v(i)=v(i)+4*x(i).*((x(i+1).^2).^(x(i).^2)).*x(i+1); v1=v; i=[(1:n)';(1:(n-1))']; j=[(1:n)';(2:n)']; s=[v0;2*v1]; H=sparse(i,j,s,n,n); H=(H+H')/2; end end