Large Sparse System of Nonlinear Equations with Jacobian

Open Live Script

This example shows how to use features of the fsolve solver to solve large sparse systems of equations effectively. The example uses the objective function, defined for a system of $n$ equations,

$\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}$

The equations to solve are $F_{i} (x) = 0$ , $1 \leq i \leq n$ . The example uses $n = 1000$ .

This objective function is simple enough that you can calculate its Jacobian analytically. As explained in Writing Vector and Matrix Objective Functions, the Jacobian $J (x)$ of a system of equations $F (x)$ is $J_{i j} (x) = \frac{\partial F_{i} (x)}{\partial x_{j}}$ . Provide this derivative as the second output of your objective function. The nlsf1 helper function at the end of this example creates the objective function $F (x)$ and its Jacobian $J (x)$ .

Solve the system of equations using the default options, which call the 'trust-region-dogleg' algorithm. Start from the point 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.

Display the solution quality and the number of function evaluations taken.

disp(norm(fval))

   2.7603e-13

disp(output.funcCount)

fsolve solves the equation accurately, but takes thousands of function evaluations to do so.

Solve the equation using the Jacobian in both the default and 'trust-region' algorithms.

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

Both algorithms take a tiny fraction of the number of function evaluations compared to the number without using the Jacobian. The default algorithm takes a few more function evaluations than the 'trust-region' algorithm, but the default algorithm reaches a more accurate answer.

See whether setting the 'PrecondBandWidth' option to 1 changes the 'trust-region' answer or efficiency. This setting causes fsolve to use a tridiagonal preconditioner, which should be effective for this tridiagonal system of equations.

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)

The 'PrecondBandWidth' option setting causes fsolve to give a slightly less accurate answer, as measured by the norm of the residual. The number of function evaluations increases slightly, from 5 to 6. The solver has fewer than ten conjugate gradient iterations as part of its solution process.

See how well fsolve performs with a diagonal preconditioner.

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)

The residual norm is slightly lower this time, and the number of function evaluations is unchanged. The number of conjugate gradient iterations increases from 8 to 19, indicating that this 'PrecondBandWidth' setting causes the solver to do more work.

Solve the equations using the 'levenberg-marquardt' algorithm.

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))

   5.2545e-15

disp(output6.funcCount)

This algorithm gives the lowest residual norm and uses only one more than the lowest number of function evaluations.

Summarize the results.

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.7603e-13     7007 
    Default+Jacobian                2.5891e-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    5.2545e-15        6

This code creates the nlsf1 function.

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

Large Sparse System of Nonlinear Equations with Jacobian

See Also

Related Topics