Main Content

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

ヤコビ スパース パターンを使った非線形方程式

解析ヤコビアンによる非線形方程式の例では、関数 nlsf1F の計算値を使ってスパースなヤコビアン J を計算します。それでは、ヤコビアンを計算するコードが使用できない場合はどうするのでしょうか?既定の設定ではユーザーが nlsf1 の中で、(options'SpecifyObjectiveGradient' オプションを true に設定して) ヤコビアンを計算しない場合、fsolvelsqnonlinlsqcurvefit は有限差分を使ってヤコビアンを近似します。

この有限差分をできるだけ効率的にするには、optionsJacobPattern をスパース行列 Jstr に設定して、ヤコビアンのスパース パターンを設定しなければなりません。すなわち、すべての x に対して、ヤコビアンの非ゼロ要素に対応する非ゼロ入力をもつスパース行列 Jstr を設定します。実際には Jstr の非ゼロ部が J の非ゼロ位置の上位集合に対応したものになりますが、一般にはスパース有限差分法の計算コストは Jstr の非ゼロ要素数と共に増加します。

スパース パターンを使うことにより、大規模問題での有限差分の計算時間を大幅に減らすことができます。スパース パターンが与えられない場合 (そしてヤコビアンが目的関数でも計算されない場合)、1,000 個の変数をもつこの問題では、有限差分を実行するコードがヤコビアンの 1,000 行 1,000 列の要素すべてを計算しようとします。しかし、この場合は 1000 行 1000 列の要素に、非ゼロ要素は 2998 個しかありません。言い換えれば、この問題を解くには、スパースのパターンが必要になります。スパースを取り扱えない場合、密行列を解こうとして、ほとんどの計算機ではメモリ不足になってしまいます。しかし、小規模な問題では、スパースの構造を与えることは必須ではありません。

事前に演算したスパース行列 Jstr をファイル nlsdat1.mat に保存しているとします。次のドライバーは、ヤコビアンなしの nlsf1 と同じ nlsf1a に適用する fsolve を呼び出します。スパース有限差分は必要に応じてスパースなヤコビ行列を評価するために計算されます。

手順 1: 目的関数値を計算するファイル nlsf1a.m を記述する

function F = nlsf1a(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;

手順 2: 方程式系を解くソルバー ルーチンを呼び出す

xstart = -ones(1000,1); 
fun = @nlsf1a;
load nlsdat1   % Get Jstr
options = optimoptions(@fsolve,'Display','iter','JacobPattern',Jstr,...
    'Algorithm','trust-region','SubproblemAlgorithm','cg');
[x,fval,exitflag,output] = fsolve(fun,xstart,options);

この場合、以下の出力が表示されます。

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          5            1011                            19
     1         10         16.1942        7.91898           2.35      
     2         15       0.0228025        1.33142          0.291      
     3         20      0.00010336      0.0433327         0.0201      
     4         25     7.37924e-07      0.0022606       0.000946      
     5         30     4.02301e-10    0.000268382       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.

またはスパースな直接的線形ソルバー (すなわちスパース QR 分解) に「complete」前提条件子を用いて選択できます。たとえば、PrecondBandWidthInf に設定した場合、スパースな直接的線形ソルバーが前処理付き共役勾配の反復の代わりに使用されます。

xstart = -ones(1000,1);
fun = @nlsf1a;
load nlsdat1   % Get Jstr
options = optimoptions(@fsolve,'Display','iter','JacobPattern',Jstr,...
 'Algorithm','trust-region','SubproblemAlgorithm','factorization');
[x,fval,exitflag,output] = fsolve(fun,xstart,options);

以下の結果が表示されます。

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          5            1011                            19
     1         10         15.9018        7.92421           1.89      
     2         15       0.0128161        1.32542         0.0746      
     3         20     1.73502e-08      0.0397923       0.000196      
     4         25     1.10732e-18    4.55495e-05       2.74e-09      

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.

スパース行列に直接適用するソルバーを使用する際には、CG の反復は存在しません。最終的な最適解および f(x) 値 (fsolve に対して f(x) は関数値の二乗和になります) は、多くの場合に PCG 法の使用時よりもゼロに近くなります。

参考

関連するトピック