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
$$f\left(x\right)=\widehat{f}\left(x\right)-\frac{1}{2}{x}^{T}V{V}^{T}x,$$
where V is a 1000-by-2 matrix. The Hessian of f is dense, but the Hessian of $$\widehat{f}$$ is sparse. If the Hessian of is $$\widehat{H}$$, then H, the Hessian of f, is
$$H=\widehat{H}-V{V}^{T}.$$
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);
Note
The function |
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 2297.63 1.41e+03 1 1084.59 6.3903 578 1 2 1084.59 100 578 3 3 1084.59 25 578 0 4 1084.59 6.25 578 0 5 1047.61 1.5625 240 0 6 761.592 3.125 62.4 2 7 761.592 6.25 62.4 4 8 746.478 1.5625 163 0 9 546.578 3.125 84.1 2 10 274.311 6.25 26.9 2 11 55.6193 11.6597 40 2 12 55.6193 25 40 3 13 22.2964 6.25 26.3 0 14 -49.516 6.25 78 1 15 -93.2772 1.5625 68 1 16 -207.204 3.125 86.5 1 17 -434.162 6.25 70.7 1 18 -681.359 6.25 43.7 2 19 -681.359 6.25 43.7 4 20 -698.041 1.5625 191 0 21 -723.959 3.125 256 7 22 -751.33 0.78125 154 3 23 -793.974 1.5625 24.4 3 24 -820.831 2.51937 6.11 3 25 -823.069 0.562132 2.87 3 26 -823.237 0.196753 0.486 3 27 -823.245 0.0621202 0.386 3 28 -823.246 0.0199951 0.11 6 29 -823.246 0.00731333 0.0404 7 30 -823.246 0.00505883 0.0185 8 31 -823.246 0.00126471 0.00268 9 32 -823.246 0.00149326 0.00521 9 33 -823.246 0.000373314 0.00091 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.3093e-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.