## Documentation Center |

The `quadprog` trust-region-reflective
method can also 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 active-set problems and 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 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.

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`HessMult`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 'HessMult' 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 HessMult option options = optimoptions(@quadprog,'Algorithm','trust-region-reflective','HessMult',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`,
and `output`. The results are

Optimization terminated: relative function value changing by less than sqrt(OPTIONS.TolFun), no negative curvature detected in current trust region model and the rate of progress (change in f(x)) is slow. fval,exitflag,output fval = -1.0538e+003 exitflag = 3 output = algorithm: 'trust-region-reflective' iterations: 18 constrviolation: 0 firstorderopt: 0.0043 cgiterations: 30 message: 'Optimization terminated: relative function value changing by less than s...'

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

In this example, `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 'HessMult' 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 HessMult option, and override the TolPCG option options = optimoptions(@quadprog,'Algorithm','trust-region-reflective',... 'HessMult',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

but the first-order optimality is further reduced.

output.firstorderopt ans = 0.0028

Was this topic helpful?