Main Content

ヤコビアンを使用した場合と使用しない場合の非線形系の解法

この例では、連立非線形方程式に導関数を与えた場合の関数評価回数の減少を示します。ベクトルと行列の目的関数の記述で説明されているように、連立方程式 F(x) のヤコビアン J(x)Jij(x)=Fi(x)xj です。この導関数を目的関数の 2 番目の出力として与えます。

たとえば、n が正の偶数値である場合、関数 multirosenbrock はローゼンブロックの関数の n 次元一般化です (制約付き非線形問題の解法、問題ベースを参照)。

F(1)=1-x1F(2)=10(x2-x12)F(3)=1-x3F(4)=10(x4-x32)F(n-1)=1-xn-1F(n)=10(xn-xn-12).

方程式系 F(x)=0 の解は、点 xi=1, i=1n です。

この目的関数では、ij の差が 1 以下の項を除く、すべてのヤコビアン項 Jij(x) が 0 です。i<n の奇数値では、非ゼロ項が次のようになります。

Jii(x)=-1J(i+1)i=-20xiJ(i+1)(i+1)=10.

この例の最後に記載されている multirosenbrock 補助関数が、目的関数 F(x) とそのヤコビアン J(x) を作成します。

i<n の奇数値の点 xi=-1.9 と、i の偶数値の点 xi=2 で始まる方程式系を解きます。n=64 と指定します。

n = 64;  
x0(1:n,1) = -1.9; 
x0(2:2:n,1) = 2;
[x,F,exitflag,output,JAC] = fsolve(@multirosenbrock,x0);
Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

計算された解 x と真の解の距離と、fsolve が解を計算するために実行した関数評価の回数を調査します。

disp(norm(x-ones(size(x))))
     0
disp(output.funcCount)
        1043

fsolve は、解を見つけるために、1000 回以上の関数評価を行っています。

今度はヤコビアンを使用して、再度、連立方程式を解きます。これを行うには、'SpecifyObjectiveGradient' オプションを true に設定します。

opts = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,F2,exitflag2,output2,JAC2] = fsolve(@multirosenbrock,x0,opts);
Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

今回も、計算された解 x2 と真の解の距離と、fsolve が解を計算するために実行した関数評価の回数を調査します。

disp(norm(x2-ones(size(x2))))
     0
disp(output2.funcCount)
    21

fsolve は、前回と同じ解を返しますが、そのための関数評価は 1,000 回以上ではなく、約 20 回です。一般的に、ヤコビアンを使用すると、関数評価の回数を減らして、ロバスト性を向上させることができますが、この例ではロバスト性の向上は見られません。

補助関数

次のコードは、補助関数 multirosenbrock を作成します。

function [F,J] = multirosenbrock(x)
% Get the problem size
n = length(x);  
if n == 0, error('Input vector, x, is empty.'); end
if mod(n,2) ~= 0
   error('Input vector, x ,must have an even number of components.');
end
% Evaluate the vector function
odds  = 1:2:n;
evens = 2:2:n;
F = zeros(n,1);
F(odds,1)  = 1-x(odds);
F(evens,1) = 10.*(x(evens)-x(odds).^2); 
% Evaluate the Jacobian matrix if nargout > 1
if nargout > 1
   c = -ones(n/2,1);    C = sparse(odds,odds,c,n,n);
   d = 10*ones(n/2,1);  D = sparse(evens,evens,d,n,n);
   e = -20.*x(odds);    E = sparse(evens,odds,e,n,n);
   J = C + D + E;
end
end

参考

関連するトピック