| Optimization Toolbox | |
Nonlinear Equations with Jacobian Sparsity Pattern
In the preceding example, the function nlsf1 computes the Jacobian J, a sparse matrix, along with the evaluation of F. What if the code to compute the Jacobian is not available? By default, if you do not indicate that the Jacobian can be computed in nlsf1 (using the Jacobian parameter in options), fsolve, lsqnonlin, and lsqcurvefit instead uses finite differencing to approximate the Jacobian.
In order for this finite differencing to be as efficient as possible, you should supply the sparsity pattern of the Jacobian, using the JacobPattern parameter in options. That is, supply a sparse matrix Jstr whose nonzero entries correspond to nonzeros of the Jacobian for all x. Indeed, the nonzeros of Jstr can correspond to a superset of the nonzero locations of J; however, in general the computational cost of the sparse finite-difference procedure will increase with the number of nonzeros of Jstr.
Providing the sparsity pattern can drastically reduce the time needed to compute the finite differencing on large problems. If the sparsity pattern is not provided (and the Jacobian is not computed in the objective function either) then, in this problem nlsfs1, the finite-differencing code attempts to compute all 1000-by-1000 entries in the Jacobian. But in this case there are only 2998 nonzeros, substantially less than the 1,000,000 possible nonzeros the finite-differencing code attempts to compute. In other words, this problem is solvable if you provide the sparsity pattern. If not, most computers run out of memory when the full dense finite-differencing is attempted. On most small problems, it is not essential to provide the sparsity structure.
Suppose the sparse matrix Jstr, computed previously, has been saved in file nlsdat1.mat. The following driver calls fsolve applied to nlsf1a, which is the same as nlsf1 except that only the function values are returned; sparse finite-differencing is used to estimate the sparse Jacobian matrix as needed.
Step 1: Write an M-file nlsf1a.m that computes the objective function values.
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;
Step 2: Call the system of equations solve routine.
xstart = -ones(1000,1); fun = @nlsf1a; load nlsdat1 % Get Jstr options = optimset('Display','iter','JacobPattern',Jstr,... 'LargeScale','on','PrecondBandWidth',1); [x,fval,exitflag,output] = fsolve(fun,xstart,options);
In this case, the output displayed is
Norm of First-order CG- Iteration Func-count f(x) step optimality Iterations 1 6 1011 1 19 0 2 11 16.0839 7.92496 1.92 1 3 16 0.0458181 1.3279 0.579 1 4 21 0.000101184 0.0631898 0.0203 2 5 26 3.16615e-007 0.00273698 0.00079 2 6 31 9.72482e-010 0.00018111 5.82e-005 2 Optimization terminated successfully: Relative function value changing by less than OPTIONS.TolFun
Alternatively, it is possible to choose a sparse direct linear solver (i.e., a sparse QR factorization) by indicating a "complete" preconditioner. I.e., if you set PrecondBandWidth to Inf, then a sparse direct linear solver is used instead of a preconditioned conjugate gradient iteration:
xstart = -ones(1000,1); fun = @nlsf1a; load nlsdat1 % Get Jstr options = optimset('Display','iter','JacobPattern',Jstr,... 'LargeScale','on','PrecondBandWidth',inf); [x,fval,exitflag,output] = fsolve(fun,xstart,options);
Norm of First-order CG- Iteration Func-count f(x) step optimality Iterations 1 6 1011 1 19 0 2 11 15.9018 7.92421 1.89 1 3 16 0.0128163 1.32542 0.0746 1 4 21 1.73538e-008 0.0397925 0.000196 1 5 26 1.13169e-018 4.55544e-005 2.76e-009 1 Optimization terminated successfully: Relative function value changing by less than OPTIONS.TolFun
When the sparse direct solvers are used, the CG iteration is 1 for that (major) iteration, as shown in the output under CG-Iterations. Notice that the final optimality and f(x) value (which for fsolve, f(x), is the sum of the squares of the function values) are closer to zero than using the PCG method, which is often the case.
| | Nonlinear Equations with Jacobian | Nonlinear Least-Squares with Full Jacobian Sparsity Pattern | |