Main Content

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

ヤコビアンを使った非線形方程式

ヤコビアンがスパースである非線形方程式系の解を求める問題を考えてみましょう。この例の問題の次元は 1000 です。ゴールは F(x) = 0 となる x を求めることです。n = 1000 を仮定して、非線形方程式は次のとおりです。

F(1)=3x12x122x2+1,F(i)=3xi2xi2xi12xi+1+1,F(n)=3xn2xn2xn1+1.

大規模な非線形方程式系 F(x) = 0 を解くには、大規模アルゴリズム (大規模アルゴリズムと中規模アルゴリズム) である fsolve で使用可能な信頼領域 Reflective 法アルゴリズムを使用してください。

手順 1: 目的関数値およびヤコビアンを計算するファイル nlsf1.m を記述する

function [F,J] = nlsf1(x)
% Evaluate the vector function
n = length(x);
F = zeros(n,1);
i = 2:(n-1);
F(i) = (3-2*x(i)).*x(i)-x(i-1)-2*x(i+1) + 1;
F(n) = (3-2*x(n)).*x(n)-x(n-1) + 1;
F(1) = (3-2*x(1)).*x(1)-2*x(2) + 1;
% Evaluate the Jacobian if nargout > 1
if nargout > 1
   d = -4*x + 3*ones(n,1); D = sparse(1:n,1:n,d,n,n);
   c = -2*ones(n-1,1); C = sparse(1:n-1,2:n,c,n,n);
   e = -ones(n-1,1); E = sparse(2:n,1:n-1,e,n,n);
   J = C + D + E;
end

手順 2: 方程式系に対するソルバー ルーチンを呼び出す

xstart = -ones(1000,1);
fun = @nlsf1; 
options = optimoptions(@fsolve,'Display','iter',...
    'Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true,'PrecondBandWidth',0);
[x,fval,exitflag,output] = fsolve(fun,xstart,options);

開始値を関数名と共に与えます。fsolve の既定方法は Trust-region-dogleg 法です。そのため信頼領域法アルゴリズムを選択するために、options 引数の 'trust-region' として 'Algorithm' を指定しなければなりません。Display オプションを 'iter' に設定すると、fsolve は反復ごとに出力を表示します。'SpecifyObjectiveGradient'true に設定すると、fsolvenlsf1.m に利用可能なヤコビ情報を使用します。

上記コマンドにより、以下が出力表示されます。

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          1            1011                            19
     1          2         16.1942        7.91898           2.35      
     2          3       0.0228027        1.33142          0.291      
     3          4     0.000103359      0.0433329         0.0201      
     4          5      7.3792e-07      0.0022606       0.000946      
     5          6     4.02299e-10    0.000268381       4.12e-05

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

線形システムの解法は、メインの反復ごとに、前処理付き共役勾配法を使用して、(近似的に ) 行われます。optionsPrecondBandWidth を 0 に設定すると、対角型をした前提条件子が使用されます (PrecondBandWidth は前処理を行う行列の帯域幅を指定します。帯域幅 0 は、行列内に対角が 1 つのみであることを意味します )。

1 次の最適性に関する値から、高速な線形収束が生じます。メインの反復に必要な共役勾配 (CG) の反復数は少なく、1000 次元の問題に対してもたかだか 5 です。このことは、例の場合、線形システムが、(収束が進行するに連れて、より多くの作業が必要であるにもかかわらず ) 必ずしも解法が困難ではないことを示唆します。

3 重主対角型の前提条件子、すなわち 3 つの対角 (または帯域幅 1) をもつ前処理行列を使用する場合、PrecondBandWidth を 値 1 に設定します。

options = optimoptions(@fsolve,'Display','iter','SpecifyObjectiveGradient',true,...
   'Algorithm','trust-region','PrecondBandWidth',1);
[x,fval,exitflag,output] = fsolve(fun,xstart,options);

この場合出力は以下になります。

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          1            1011                            19
     1          2         16.0839        7.92496           1.92      
     2          3       0.0458181         1.3279          0.579      
     3          4     0.000101184      0.0631898         0.0203      
     4          5     3.16615e-07     0.00273698        0.00079      
     5          6     9.72481e-10     0.00018111       5.82e-05

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

同数の反復が生じているにもかかわらず、PCG の反復数が少なくなっていることに注意してください。これは、一回の反復で実行される作業が少なくなっていることを示します。詳細は、前処理付き共役勾配法を参照してください。

PrecondBandWidthInf に設定する (既定の設定) と、ソルバーは PCG ではなくコレスキー分解を使用します。

参考

関連するトピック