## Documentation Center |

The `fmincon` `interior-point` and `trust-region-reflective` algorithms,
and the `fminunc` `trust-region` algorithm
can solve problems where the Hessian is dense but structured. For
these problems, `fmincon` and `fminunc` do
not compute *H*Y* with the Hessian *H* directly,
because forming *H* would be memory-intensive. Instead,
you must provide `fmincon` or `fminunc` with
a function that, given a matrix *Y* and information
about *H*, computes *W* = *H*Y*.

In this example, the objective function is nonlinear and linear
equalities exist so `fmincon` is used. The description
applies to the trust-region reflective algorithm; the `fminunc` `trust-region` algorithm
is similar. For the interior-point algorithm, see the `'HessMult'` option
in Hessian. The objective function
has the structure

where *V* is a 1000-by-2 matrix. The Hessian
of *f* is dense, but the Hessian of
is sparse. If the Hessian of
is
, then *H*,
the Hessian of *f*, is

To avoid excessive memory usage that could happen by working
with *H* directly, the example provides a Hessian
multiply function, `hmfleq1`. This function, when
passed a matrix `Y`, uses sparse matrices `Hinfo`,
which corresponds to
, and `V` to compute the Hessian
matrix product

W = H*Y = (Hinfo - V*V')*Y

In this example, the Hessian multiply function needs
and `V` to
compute the Hessian matrix product. `V` is a constant,
so you can capture `V` in a function handle to an
anonymous function.

However,
is not a constant and must be computed at the
current `x`. You can do this by computing
in the objective
function and returning
as `Hinfo` in the third output
argument. By using `optimoptions` to set the `'Hessian'` options
to `'on'`, `fmincon` knows to get
the `Hinfo` value from the objective function and
pass it to the Hessian multiply function `hmfleq1`.

The example passes `brownvv` to `fmincon` as
the objective function. The `brownvv.m``brownvv.m` file
is long and is not included here. You can view the code with the command

type brownvv

Because `brownvv` computes the gradient and
part of the Hessian as well as the objective function, the example
(Step
3) uses `optimoptions` to
set the `GradObj` and `Hessian` options
to `'on'`.

Now, define a function `hmfleq1` that uses `Hinfo`,
which is computed in `brownvv`, and `V`,
which you can capture in a function handle to an anonymous function,
to compute the Hessian matrix product `W` where `W = H*Y = (Hinfo - V*V')*Y`. This function must
have the form

W = hmfleq1(Hinfo,Y)

The first argument must be the same as the third argument returned
by the objective function `brownvv`. The second argument
to the Hessian multiply function is the matrix `Y` (of `W
= H*Y`).

Because `fmincon` expects the second argument `Y` to
be used to form the Hessian matrix product, `Y` is
always a matrix with `n` rows where `n` is
the number of dimensions in the problem. The number of columns in `Y` can
vary. Finally, you can use a function handle to an anonymous function
to capture V, so V can be the third argument to `'hmfleqq'`.

function W = hmfleq1(Hinfo,Y,V); %HMFLEQ1 Hessian-matrix product function for BROWNVV objective. % W = hmfleq1(Hinfo,Y,V) computes W = (Hinfo-V*V')*Y % where Hinfo is a sparse matrix computed by BROWNVV % and V is a 2 column matrix. W = Hinfo*Y - V*(V'*Y);

Load the problem parameter, `V`, and the sparse
equality constraint matrices, `Aeq` and `beq`,
from `fleq1.mat`, which is available in the `optimdemos` folder.
Use `optimoptions` to set the `GradObj` and `Hessian` options
to `'on'` and to set the `HessMult` option
to a function handle that points to `hmfleq1`. Call `fmincon` with
objective function `brownvv` and with `V` as
an additional parameter:

function [fval, exitflag, output, x] = runfleq1 % RUNFLEQ1 demonstrates 'HessMult' option for FMINCON with linear % equalities. problem = load('fleq1'); % Get V, Aeq, beq V = problem.V; Aeq = problem.Aeq; beq = problem.beq; n = 1000; % problem dimension xstart = -ones(n,1); xstart(2:2:n,1) = ones(length(2:2:n),1); % starting point options = optimoptions(@fmincon,'Algorithm','trust-region-reflective','GradObj','on', ... 'Hessian','user-supplied','HessMult',@(Hinfo,Y)hmfleq1(Hinfo,Y,V),'Display','iter', ... 'TolFun',1e-9); [x,fval,exitflag,output] = fmincon(@(x)brownvv(x,V),xstart,[],[],Aeq,beq,[],[], ... [],options);

To run the preceding code, enter

[fval,exitflag,output,x] = runfleq1;

Because the iterative display was set using `optimoptions`,
this command generates the following iterative display:

Norm of First-order Iteration f(x) step optimality CG-iterations 0 1997.07 916 1 1072.57 6.31716 465 1 2 480.247 8.19711 201 2 3 136.982 10.3039 78.1 2 4 44.416 9.04685 16.7 2 5 44.416 100 16.7 2 6 44.416 25 16.7 0 7 -9.05631 6.25 52.9 0 8 -317.437 12.5 91.7 1 9 -405.381 12.5 1.11e+003 1 10 -451.161 3.125 327 4 11 -482.688 0.78125 303 5 12 -547.427 1.5625 187 5 13 -610.42 1.5625 251 7 14 -711.522 1.5625 143 3 15 -802.98 3.125 165 3 16 -820.431 1.13329 32.9 3 17 -822.996 0.492813 7.61 2 18 -823.236 0.223154 1.68 3 19 -823.245 0.056205 0.529 3 20 -823.246 0.0150139 0.0342 5 21 -823.246 0.00479085 0.0152 7 22 -823.246 0.00353697 0.00828 9 23 -823.246 0.000884242 0.005 9 24 -823.246 0.0012715 0.00125 9 25 -823.246 0.000317876 0.0025 9 Local minimum possible. fmincon stopped because the final change in function value relative to its initial value is less than the selected value of the function tolerance.

Convergence is rapid for a problem of this size with the PCG iteration cost increasing modestly as the optimization progresses. Feasibility of the equality constraints is maintained at the solution.

problem = load('fleq1'); % Get V, Aeq, beq V = problem.V; Aeq = problem.Aeq; beq = problem.beq; norm(Aeq*x-beq,inf) ans = 2.4869e-14

In this example, `fmincon` cannot use `H` to
compute a preconditioner because `H` only exists
implicitly. Instead of `H`, `fmincon` uses `Hinfo`,
the third argument returned by `brownvv`, to compute
a preconditioner. `Hinfo` is a good choice because
it is the same size as `H` and approximates `H` to
some degree. If `Hinfo` were not the same size as `H`, `fmincon` would
compute a preconditioner based on some diagonal scaling matrices determined
from the algorithm. Typically, this would not perform as well.

Was this topic helpful?