Step 1: Decide what part of H to pass to quadprog as the first argument.
Step 2: Write a function to compute Hessian-matrix products for H.
Step 3: Call a quadratic minimization routine with a starting point.
The quadprog
trust-region-reflective
method can solve large problems where the Hessian is dense but structured.
For these problems, quadprog
does not compute H*Y with
the Hessian H directly, as it does for trust-region-reflective
problems with sparse H, because forming H would
be memory-intensive. Instead, you must provide quadprog
with
a function that, given a matrix Y and information
about H, computes W = H*Y.
In this example, the Hessian matrix H
has
the structure H = B + A*A'
where B
is
a sparse 512-by-512 symmetric matrix, and A
is
a 512-by-10 sparse matrix composed of a number of dense columns. To
avoid excessive memory usage that could happen by working with H
directly
because H
is dense, the example provides a Hessian
multiply function, qpbox4mult
. This function, when
passed a matrix Y
, uses sparse matrices A
and B
to
compute the Hessian matrix product W = H*Y = (B + A*A')*Y
.
In the first part of this example, the matrices A
and B
need
to be provided to the Hessian multiply function qpbox4mult
.
You can pass one matrix as the first argument to quadprog
,
which is passed to the Hessian multiply function. You can use a nested
function to provide the value of the second matrix.
The second part of the example shows how to tighten the TolPCG
tolerance
to compensate for an approximate preconditioner instead of an exact H
matrix.
Either A
or B
can be passed
as the first argument to quadprog
. The example
chooses to pass B
as the first argument because
this results in a better preconditioner (see Preconditioning).
quadprog(B,f,[],[],[],[],l,u,xstart,options)
Now, define a function runqpbox4
that
Contains a nested function qpbox4mult
that
uses A
and B
to compute the
Hessian matrix product W
, where W = H*Y = (B + A*A')*Y
. The nested function must have the form
W = qpbox4mult(Hinfo,Y,...)
The first two arguments Hinfo
and Y
are
required.
Loads the problem parameters from qpbox4.mat
.
Uses optimoptions
to
set the HessianMultiplyFcn
option to a function
handle that points to qpbox4mult
.
Calls quadprog
with B
as
the first argument.
The first argument to the nested function qpbox4mult
must
be the same as the first argument passed to quadprog
,
which in this case is the matrix B.
The second argument to qpbox4mult
is the
matrix Y
(of W = H*Y
).
Because quadprog
expects 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. The function qpbox4mult
is nested so that
the value of the matrix A
comes from the outer
function. Optimization Toolbox™ software includes the runqpbox4.m
file.
function [fval, exitflag, output, x] = runqpbox4 %RUNQPBOX4 demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds. problem = load('qpbox4'); % Get xstart, u, l, B, A, f xstart = problem.xstart; u = problem.u; l = problem.l; B = problem.B; A = problem.A; f = problem.f; mtxmpy = @qpbox4mult; % function handle to qpbox4mult nested function % Choose algorithm and the HessianMultiplyFcn option options = optimoptions(@quadprog,'Algorithm','trust-region-reflective','HessianMultiplyFcn',mtxmpy); % Pass B to qpbox4mult via the H argument. Also, B will be used in % computing a preconditioner for PCG. [x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options); function W = qpbox4mult(B,Y) %QPBOX4MULT Hessian matrix product with dense structured Hessian. % W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where % INPUT: % B - sparse square matrix (512 by 512) % Y - vector (or matrix) to be multiplied by B + A'*A. % VARIABLES from outer function runqpbox4: % A - sparse matrix with 512 rows and 10 columns. % % OUTPUT: % W - The product (B + A*A')*Y. % % Order multiplies to avoid forming A*A', % which is large and dense W = B*Y + A*(A'*Y); end end
To call the quadratic minimizing routine contained in runqpbox4
,
enter
[fval,exitflag,output] = runqpbox4;
to run the preceding code. Then display the values for fval
, exitflag
, output.iterations
,
and output.cgiterations
.
fval,exitflag,output.iterations, output.cgiterations fval = -1.0538e+03 exitflag = 3 ans = 18 ans = 30
After 18 iterations with a total of 30 PCG iterations, the function value is reduced to
fval fval = -1.0538e+003
and the first-order optimality is
output.firstorderopt ans = 0.0043
Sometimes quadprog
cannot use H
to
compute a preconditioner because H
only exists
implicitly. Instead, quadprog
uses B
,
the argument passed in instead of H
, to compute
a preconditioner. B
is a good choice because it
is the same size as H
and approximates H
to
some degree. If B
were not the same size as H
, quadprog
would
compute a preconditioner based on some diagonal scaling matrices determined
from the algorithm. Typically, this would not perform as well.
Because the preconditioner is more approximate than when H
is
available explicitly, adjusting the TolPCG
parameter
to a somewhat smaller value might be required. This example is the
same as the previous one, but reduces TolPCG
from
the default 0.1 to 0.01.
function [fval, exitflag, output, x] = runqpbox4prec %RUNQPBOX4PREC demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds. problem = load('qpbox4'); % Get xstart, u, l, B, A, f xstart = problem.xstart; u = problem.u; l = problem.l; B = problem.B; A = problem.A; f = problem.f; mtxmpy = @qpbox4mult; % function handle to qpbox4mult nested function % Choose algorithm, the HessianMultiplyFcn option, and override the TolPCG option options = optimoptions(@quadprog,'Algorithm','trust-region-reflective',... 'HessianMultiplyFcn',mtxmpy,'TolPCG',0.01); % Pass B to qpbox4mult via the H argument. Also, B will be used in % computing a preconditioner for PCG. % A is passed as an additional argument after 'options' [x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options); function W = qpbox4mult(B,Y) %QPBOX4MULT Hessian matrix product with dense structured Hessian. % W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where % INPUT: % B - sparse square matrix (512 by 512) % Y - vector (or matrix) to be multiplied by B + A'*A. % VARIABLES from outer function runqpbox4prec: % A - sparse matrix with 512 rows and 10 columns. % % OUTPUT: % W - The product (B + A*A')*Y. % % Order multiplies to avoid forming A*A', % which is large and dense W = B*Y + A*(A'*Y); end end
Now, enter
[fval,exitflag,output] = runqpbox4prec;
to run the preceding code. After 18 iterations and 50 PCG iterations, the function value has the same value to five significant digits
fval fval = -1.0538e+003
and the first-order optimality is essentially the same.
output.firstorderopt ans = 0.0043
Note
Decreasing |