%% Minimization with Gradient and Hessian Sparsity Pattern
% This example shows how to solve a nonlinear minimization problem with a tridiagonal Hessian matrix approximated by sparse finite differences instead of explicit computation.
% The problem is to find x to minimize
% \[f(x) = \sum\limits_{i = 1}^{n - 1} {\left( {{{\left( {x_i^2} \right)}^{\left( {x_{i + 1}^2 + 1} \right)}} + {{\left( {x_{i + 1}^2} \right)}^{\left( {x_i^2 + 1} \right)}}} \right)} ,\]
% where n = 1000.

n = 1000;

% To use the trust-region method in fminunc, you must compute the gradient in the objective function; it is not optional as in the quasi-newton method.
% The helper function brownfg at the end of this example computes the objective function and gradient.
% To allow efficient computation of the sparse finite-difference approximation of the Hessian matrix H(x), the sparsity structure of H must be predetermined. In this case, the structure Hstr, a sparse matrix, is available in the file brownhstr.mat. Using the spy command, you can see that Hstr is, indeed, sparse (only 2998 nonzeros).

load brownhstr
spy(Hstr)

%
% Set the HessPattern option to Hstr using optimoptions. When such a large problem has obvious sparsity structure, not setting the HessPattern option uses a great amount of memory and computation unnecessarily, because fminunc attempts to use finite differencing on a full Hessian matrix of one million nonzero entries.
% To use the Hessian sparsity pattern, you must use the trust-region algorithm of fminunc. This algorithm also requires you to set the SpecifyObjectiveGradient option to true using optimoptions.

options = optimoptions(@fminunc,'Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true,'HessPattern',Hstr);

%
% Set the objective function to @brownfg. Set the initial point to –1 for odd x components and +1 for even x components.

xstart = -ones(n,1); 
xstart(2:2:n,1) = 1;
fun = @brownfg;

% Solve the problem by calling fminunc using the initial point xstart and options options.

[x,fval,exitflag,output] = fminunc(fun,xstart,options);

%
% Examine the solution and solution process.

disp(fval)
disp(exitflag)
disp(output)

% The function f(x) is a sum of powers of squares and, therefore, is nonnegative. The solution fval is nearly zero, so it is clearly a minimum. The exit flag 1 also indicates that fminunc finds a solution. The output structure shows that fminunc takes only seven iterations to reach the solution.
% Display the largest and smallest elements of the solution.

disp(max(x))
disp(min(x))

% The solution is near the point where all elements of x = 0.
%
%% Helper Function
% This code creates the brownfg helper function.

function [f,g] = brownfg(x)
% BROWNFG Nonlinear minimization test problem
% 
% Evaluate the function
n=length(x); y=zeros(n,1);
i=1:(n-1);
y(i)=(x(i).^2).^(x(i+1).^2+1) + ...
        (x(i+1).^2).^(x(i).^2+1);
  f=sum(y);
% Evaluate the gradient if nargout > 1
  if nargout > 1
     i=1:(n-1); g = zeros(n,1);
     g(i) = 2*(x(i+1).^2+1).*x(i).* ...
              ((x(i).^2).^(x(i+1).^2))+ ...
              2*x(i).*((x(i+1).^2).^(x(i).^2+1)).* ...
              log(x(i+1).^2);
     g(i+1) = g(i+1) + ...
              2*x(i+1).*((x(i).^2).^(x(i+1).^2+1)).* ...
              log(x(i).^2) + ...
              2*(x(i).^2+1).*x(i+1).* ...
              ((x(i+1).^2).^(x(i).^2));
  end
end

% Copyright 2019 The MathWorks, Inc.