Quadratic Minimization with Dense, Structured Hessian
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.
Take advantage of a structured Hessian
Thequadprog
trust-region-reflective method can solve large problems where the Hessian is dense but structured. For these problems,quadprog
does not computeH*Ywith the HessianHdirectly, as it does for trust-region-reflective problems with sparseH, because formingHwould be memory-intensive. Instead, you must providequadprog
函数,给定一个矩阵Yand information aboutH, computesW=H*Y.
In this example, the Hessian matrixH
has the structureH = B + A*A'
whereB
is a sparse 512-by-512 symmetric matrix, andA
is a 512-by-10 sparse matrix composed of a number of dense columns. To avoid excessive memory usage that could happen by working withH
directly becauseH
is dense, the example provides a Hessian multiply function,qpbox4mult
. This function, when passed a matrixY
, uses sparse matricesA
andB
to compute the Hessian matrix productW = H*Y = (B + A*A')*Y
.
In the first part of this example, the matricesA
andB
need to be provided to the Hessian multiply functionqpbox4mult
. You can pass one matrix as the first argument toquadprog
, 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 theTolPCG
tolerance to compensate for an approximate preconditioner instead of an exactH
matrix.
Step 1: Decide what part of H to pass to quadprog as the first argument.
EitherA
orB
can be passed as the first argument toquadprog
. The example chooses to passB
as the first argument because this results in a better preconditioner (seePreconditioning).
quadprog(B,f,[],[],[],[],l,u,xstart,options)
Step 2: Write a function to compute Hessian-matrix products for H.
Now, define a functionrunqpbox4
that
Contains a nested function
qpbox4mult
that usesA
andB
to compute the Hessian matrix productW
, whereW = H*Y = (B + A*A')*Y
. The nested function must have the formW = qpbox4mult(Hinfo,Y,...)
The first two arguments
Hinfo
andY
are required.Loads the problem parameters from
qpbox4.mat
.Uses
optimoptions
to set theHessianMultiplyFcn
option to a function handle that points toqpbox4mult
.Calls
quadprog
withB
as the first argument.
The first argument to the nested functionqpbox4mult
must be the same as the first argument passed toquadprog
, which in this case is the matrix B.
The second argument toqpbox4mult
is the matrixY
(ofW = H*Y
). Becausequadprog
expectsY
to be used to form the Hessian matrix product,Y
is always a matrix withn
rows, wheren
is the number of dimensions in the problem. The number of columns inY
can vary. The functionqpbox4mult
is nested so that the value of the matrixA
comes from the outer function. Optimization Toolbox™ software includes therunqpbox4.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
Step 3: Call a quadratic minimization routine with a starting point.
To call the quadratic minimizing routine contained inrunqpbox4
, enter
[fval,exitflag,output] = runqpbox4;
to run the preceding code. Then display the values forfval
,exitflag
,output.iterations
, andoutput.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
Preconditioning
Sometimesquadprog
cannot useH
to compute a preconditioner becauseH
only exists implicitly. Instead,quadprog
usesB
, the argument passed in instead ofH
, to compute a preconditioner.B
is a good choice because it is the same size asH
and approximatesH
to some degree. IfB
were not the same size asH
,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 whenH
is available explicitly, adjusting theTolPCG
parameter to a somewhat smaller value might be required. This example is the same as the previous one, but reducesTolPCG
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, fxstart = 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 optionoptions = 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);functionW = 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 denseW = B*Y + A*(A'*Y);endend
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
DecreasingTolPCG
too much can substantially increase the number of PCG iterations.