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 formW = 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

Decreasing `TolPCG`

too much can substantially
increase the number of PCG iterations.