Main Content

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\le i\le 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}_{ij}(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.8577e-13

disp(output.funcCount)

7007

The equation is solved accurately, but `fsolve`

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)

6

disp(output4.cgiterations)

8

This option setting causes `fsolve`

to give a slightly less accurate answer as measured by the norm of the residual. The number of function evaluations increased slightly, from 5 to 6. The solver had 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)

6

disp(output5.cgiterations)

19

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

disp(norm(fval6))

7.5854e-15

disp(output6.funcCount)

6

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.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 7.5854e-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