Documentation

The goal in this problem is to minimize the nonlinear function

$$f(x)=1+{\displaystyle \sum _{i=1}^n{\left|\left(3-2x_i\right)x_i-x_{i-1}-x_{i+1}+1\right|}^p}+{\displaystyle \sum _{i=1}^{n/2}{\left|x_i+x_{i+n/2}\right|}^p},$$

such that -10.0 ≤ x_i ≤ 10.0, where n is 800 (n should be a multiple of 4), p = 7/3, and x₀ = x_n + 1 = 0.

Step 1: Write a file tbroyfg.m that computes the objective function and the gradient of the objective

The tbroyfg.m file computes the function value and gradient. This file is long and is not included here. You can see the code for this function using the command

type tbroyfg

The sparsity pattern of the Hessian matrix has been predetermined and stored in the file tbroyhstr.mat. The sparsity structure for the Hessian of this problem is banded, as you can see in the following spy plot.

load tbroyhstr
spy(Hstr)

In this plot, the center stripe is itself a five-banded matrix. The following plot shows the matrix more clearly:

spy(Hstr(1:20,1:20))

Use optimoptions to set the HessPattern parameter to Hstr. When a problem as large as this has obvious sparsity structure, not setting the HessPattern parameter requires a huge amount of unnecessary memory and computation. This is because fmincon attempts to use finite differencing on a full Hessian matrix of 640,000 nonzero entries.

You must also set the GradObj parameter to 'on' using optimoptions, since the gradient is computed in tbroyfg.m. Then execute fmincon as shown in Step 2.

Step 2: Call a nonlinear minimization routine with a starting point xstart.

fun = @tbroyfg;
load tbroyhstr          % Get Hstr, structure of the Hessian
n = 800;
xstart = -ones(n,1); xstart(2:2:n) = 1;
lb = -10*ones(n,1); ub = -lb;
options = optimoptions('fmincon','GradObj','on','HessPattern',Hstr,...
    'Algorithm','trust-region-reflective'); 

[x,fval,exitflag,output] = ... 
   fmincon(fun,xstart,[],[],[],[],lb,ub,[],options);

After seven iterations, the exitflag, fval, and output values are

exitflag =
     3

fval =
  270.4790

output = 
         iterations: 7
          funcCount: 8
           stepsize: 8.2073e-04
       cgiterations: 18
      firstorderopt: 0.0163
          algorithm: 'trust-region-reflective'
            message: 'Local minimum possible.…'
    constrviolation: 0

For bound constrained problems, the first-order optimality is the infinity norm of v.*g, where v is defined as in Box Constraints, and g is the gradient.

Because of the five-banded center stripe, you can improve the solution by using a five-banded preconditioner instead of the default diagonal preconditioner. Using the optimoptions function, reset the PrecondBandWidth parameter to 2 and solve the problem again. (The bandwidth is the number of upper (or lower) diagonals, not counting the main diagonal.)

fun = @tbroyfg;
load tbroyhstr          % Get Hstr, structure of the Hessian
n = 800;
xstart = -ones(n,1); xstart(2:2:n,1) = 1;
lb = -10*ones(n,1); ub = -lb;
options = optimoptions('fmincon','GradObj','on','HessPattern',Hstr, ...
    'Algorithm','trust-region-reflective','PrecondBandWidth',2); 
[x,fval,exitflag,output] = ...
   fmincon(fun,xstart,[],[],[],[],lb,ub,[],options); 

The number of iterations actually goes up by two; however the total number of CG iterations drops from 18 to 15. The first-order optimality measure is reduced by a factor of 1e-3:

exitflag =
     3

fval =
  270.4790

output = 
         iterations: 9
          funcCount: 10
           stepsize: 2.4512e-05
       cgiterations: 15
      firstorderopt: 7.5340e-05
          algorithm: 'trust-region-reflective'
            message: 'Local minimum possible.…'
    constrviolation: 0