ヤコビアンを使用した非線形方程式の大規模スパース系

この例では、fsolve ソルバーの機能を使用して、方程式の大規模スパース系を効率的に解く方法を示します。この例では、 $n$ 個の方程式の系として定義された目的関数を使用します。

$\begin{array}{l} F (1) = 3 x_{1} - 2 x_{1}^{2} - 2 x_{2} + 1, \\ F (i) = 3 x_{i} - 2 x_{i}^{2} - x_{i - 1} - 2 x_{i + 1} + 1, \\ F (n) = 3 x_{n} - 2 x_{n}^{2} - x_{n - 1} + 1 . \end{array}$

解く方程式は、 $F_{i} (x) = 0$ , $1 \leq i \leq n$ です。この例では、 $n = 1000$ を使用します。

この目的関数は、そのヤコビアンを分析的に計算できるほど十分単純です。ベクトルと行列の目的関数の記述で説明されているように、連立方程式 $F (x)$ のヤコビアン $J (x)$ は $J_{i j} (x) = \frac{\partial F_{i} (x)}{\partial x_{j}}$ です。この導関数を目的関数の 2 番目の出力として与えます。この例の最後に記載されている nlsf1 補助関数が、目的関数 $F (x)$ とそのヤコビアン $J (x)$ を作成します。

'trust-region-dogleg' アルゴリズムを呼び出す既定のオプションを使用して、連立方程式を解きます。点 xstart(i) = -1 から始めます。

n = 1000;
xstart = -ones(n,1);
fun = @nlsf1; 
[x,fval,exitflag,output] = fsolve(fun,xstart);

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.

解の質と実行された関数評価の回数を表示します。

disp(norm(fval))

   2.8577e-13

disp(output.funcCount)

fsolve は方程式を正確に解いていますが、そのために何千回もの関数評価を行っています。

既定のアルゴリズムと 'trust-region' アルゴリズムの両方でヤコビアンを使用して方程式を解きます。

options = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,fval2,exitflag2,output2] = fsolve(fun,xstart,options);

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.

options.Algorithm = 'trust-region';
[x3,fval3,exitflag3,output3] = fsolve(fun,xstart,options);

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.

disp([norm(fval2),norm(fval3)])

   1.0e-08 *

    0.0000    0.1065

disp([output2.funcCount,output3.funcCount])

     7     5

両方のアルゴリズムで、ヤコビアンを使用しなかった場合の回数と比べて、関数評価の回数が大幅に減っています。既定のアルゴリズムの方が 'trust-region' アルゴリズムよりも数回多く関数評価を行っていますが、既定のアルゴリズムの方がより正確な答えを出しています。

'PrecondBandWidth' オプションを 1 に設定した場合、'trust-region' の答えが変わるかあるいは効率的になるかを見てみます。この設定によって、fsolve は三重前提条件子を使用するようになります。これは、この三重対角連立方程式には有効なはずです。

options.PrecondBandWidth = 1;
[x4,fval4,exitflag4,output4] = fsolve(fun,xstart,options);

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.

disp(norm(fval4))

   3.1185e-05

disp(output4.funcCount)

disp(output4.cgiterations)

'PrecondBandWidth' オプションの設定によって、fsolve は、残差のノルムで測定した場合よりもわずかに不正確な答えを返します。関数評価の回数は、5 回から 6 回へとわずかに増えます。ソルバーの解法プロセスの一部としての共役勾配反復は 10 回未満です。

対角前提条件子を使用するとどのくらい fsolve のパフォーマンスが改善するかを確認します。

options.PrecondBandWidth = 0;
[x5,fval5,exitflag5,output5] = fsolve(fun,xstart,options);

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.

disp(norm(fval5))

   2.0057e-05

disp(output5.funcCount)

disp(output5.cgiterations)

今回は、残差ノルムがわずかに減少し、関数評価の回数は変化しません。共役勾配反復の回数は 8 回から 19 回に増加しますが、これは、この 'PrecondBandWidth' 設定によってソルバーの作業量が増えたことを示しています。

'levenberg-marquardt' アルゴリズムを使用して、方程式を解きます。

options = optimoptions('fsolve','SpecifyObjectiveGradient',true,'Algorithm','levenberg-marquardt');
[x6,fval6,exitflag6,output6] = fsolve(fun,xstart,options);

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.

disp(norm(fval6))

   8.0029e-15

disp(output6.funcCount)

このアルゴリズムでは、残差ノルムが最も少なくなり、最も少ない回数よりも 1 回多いだけの関数評価を使用します。

結果をまとめます。

t = table([norm(fval);norm(fval2);norm(fval3);norm(fval4);norm(fval5);norm(fval6)],...
    [output.funcCount;output2.funcCount;output3.funcCount;output4.funcCount;output5.funcCount;output6.funcCount],...
    'VariableNames',["Residual" "Fevals"],...
    'RowNames',["Default" "Default+Jacobian" "Trust-Region+Jacobian" "Trust-Region+Jacobian,BW=1" "Trust-Region+Jacobian,BW=0" "Levenberg-Marquardt+Jacobian"])

t=6×2 table
                                     Residual     Fevals
                                    __________    ______

    Default                         2.8577e-13     7007 
    Default+Jacobian                2.5886e-13        7 
    Trust-Region+Jacobian           1.0646e-09        5 
    Trust-Region+Jacobian,BW=1      3.1185e-05        6 
    Trust-Region+Jacobian,BW=0      2.0057e-05        6 
    Levenberg-Marquardt+Jacobian    8.0029e-15        6

このコードは関数 nlsf1 を作成します。

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
end

参考

fsolve

ヤコビアンを使用した非線形方程式の大規模スパース系

参考

関連するトピック