Ask the Expert: Sparse Nonlinear Systems using e04vh

"I need to solve a large sparse nonlinear system. How can NAG help me?"

There are currently no NAG routines available specifically for solving large sparse systems of nonlinear equations. However, it is relatively straightforward to adapt the general sparse nonlinear optimizer from the "Minimizing or Maximizing a Function" (E04) Chapter to solve sparse nonlinear systems; here the solution of such systems amounts to finding a feasible point for an optimization problem whose set of nonlinear constraints is the nonlinear system to be solved. (Functionality for dense systems exists in the "Roots of One or More Transcendental Equations" (C05) Chapter. For example, in the NAG Toolbox for MATLAB Mark 21, the relevant functions are the c05n* and the c05p* families: see C05 Chapter Contents. These dense solvers are unsuitable for systems with small percentages of nonzero Jacobian elements.)

The demonstration below shows how to solve a large sparse nonlinear system using the function e04vh from the NAG Toolbox for MATLAB Mark 21.

Consider the following problem, sometimes called "Broyden's Tridiagonal system":

for a given system size n, find x=(x1,...,xn) such that

Fi(x) = (3 - 2xi)xi - xi-1 - 2xi+1 + 1 = 0, for i = 1,...,n, defining x0 = xn+1 = 0

using the initial guess x = (-1,...,-1).

Note that there are 3*n - 2 nonzero elements in the Jacobian of F, from a total of n² elements.

Before solving this system using the NAG Toolbox for MATLAB we need to convert the mathematical formulation of the system into the special format required by e04vh (this is not compulsory but is required for optimal performance). How to do this is covered in more detail in the e04vh document, but essentially we split the system into those linear variables that appear purely linearly and those that appear at all nonlinearly:

Fi(x) = fi(x) + ai⋅x

from which we have

fi(x) = (3 - 2xi)xi

ai = -ei-1 - 2ei+1

where e_j is the jth unit vector in n-space.

Now we're ready to solve the system. The very first step when using e04vh is to initialize the solver's communication arrays

% First, call e04vg to initialize e04vh
[cw, iw, rw, ifail] = e04vg();

Then we can set the problem name and size, and also the sparse matrix whose columns are the vectors a_i. The problem name will appear in monitoring output if monitoring is enabled.

% Problem name in printout: Broyden's Tridigonal system
prob = 'BroyTrid';

% Problem size
nf = 1000;

% Set up tridiagonal A and 'sparsify'
e = ones(nf, 1);
a = spdiags([-e -2*e], [-1, 1], double(nf), double(nf));
[iafun, javar, a] = find(a);
iafun = int32(iafun);
javar = int32(javar);
nea = int32(length(a));

MATLAB's find allows us easily to obtain the coordinate arrays iafun and javar and the packed representation of the dense matrix in the form required by e04vh.

Next, set up the sparsity pattern for the Jacobian of the nonlinear part, f, of F

% Set up the sparsity pattern of the nonlinear Jacobian: diagonal
neg = nf;
igfun = int32((1:neg));
jgvar = igfun;

e04vj can be used to determine the sparsity pattern if it is not known.

The input variable n in e04vh sets the number of variables in the problem: in this case this equals nf. We also need to set objadd:

n = nf;
objadd = 0;

The bounds on the constraints in e04vh (flow and fupp) come from the specification of the system F, after rearranging so that constant terms appear on the righthand side of the definition of each F_i. There are no bounds on the variables x_i.

% The initial point for the optimizer
x = -ones(n, 1);

% The lower and upper bounds for the equations
flow = -ones(nf, 1);
fupp = flow;

% The variables are unbounded
infbnd = 1e20;
xlow = -infbnd*ones(n, 1);
xupp = -xlow;

We can flag this particular optimization problem as a feasible-point problem using objrow. We're also going to be using a cold start, so most of the input data intended for use during warm starts can just be set to zero (but must still be defined to have the correct sizes and types)

% Feasible-point problem
objrow = int32(0);

% Cold start
start = int32(0);
xstate = zeros(n, 1, 'int32');
f = zeros(nf, 1);
fstate = zeros(nf, 1, 'int32');
fmul = zeros(nf, 1);
ns = int32(0);

Arrays xnames and fnames can be used to label the variables and constraints in e04vh's monitoring output, if enabled.

% No variable names in printout
xnames = {'None'};
fnames = {'None'};

Finally, we can call the solver

[xOut, xstateOut, xmul, fOut, fstateOut, fmulOut, nsOut, ninf, sinf, ...
    cwOut, iwOut, rwOut, cuser, ifail] = ...
e04vh(start, objadd, objrow, prob, 'e04vh_usrfun', iafun, javar, ...
      nea, a, igfun, jgvar, neg, xlow, xupp, xnames, flow, fupp, ...
      fnames, x, xstate, f, fstate, fmul, ns, cw, iw, rw);

where our usrfun e04vh_usrfun looks like

function [status, f, g, user] = ...
    e04vh_usrfun(status, n, x, needf, nf, f, needg, leng, g, user)

  if (needf > 0)
      f(1:nf) = (3 - 2*x(1:nf)).*x(1:nf);
  end

  if (needg > 0)
      g(1:nf) = 3 - 4*x(1:nf);
  end

Expanded versions of this MATLAB source can be downloaded from http://www.nag.co.uk/nagnews/e04vh.zip

---

For specific technical advice in using NAG's products, please contact our technical experts.

Return to Technical Tips & Hints index page.