# Quadratic Minimization with Dense, Structured Hessian

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.

### Step 1: Decide what part of H to pass to `quadprog` as the first argument.

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).

### Step 2: Write a function to compute Hessian-matrix products for H.

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 following form, where the first two arguments `Hinfo` and `Y` are required:

```W = qpbox4mult(Hinfo,Y,...) ```
• 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. The `runqpbox4.m` file appears at the end of this example.

### Step 3: Call a quadratic minimization routine with a starting point.

To call the quadratic minimizing routine contained in `runqpbox4`, enter

`[fval,exitflag,output] = runqpbox4;`
```Local minimum possible. quadprog stopped because the relative change in function value is less than the sqrt of the function tolerance, the rate of change in the function value is slow, and no negative curvature was detected. ```

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+03 ```

and the first-order optimality is

`output.firstorderopt`
```ans = 0.0043 ```

### Preconditioning

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. The `runqpbox4prec` function is listed at the end of this example.

`[fval,exitflag,output] = runqpbox4prec;`
```Local minimum possible. quadprog stopped because the relative change in function value is less than the sqrt of the function tolerance, the rate of change in the function value is slow, and no negative curvature was detected. ```

After 18 iterations and 50 PCG iterations, the function value has the same value to five significant digits

`fval`
```fval = -1.0538e+03 ```

and the first-order optimality is essentially the same:

`output.firstorderopt`
```ans = 0.0028 ```

Note: Decreasing `TolPCG` too much can substantially increase the number of PCG iterations.

### Helper Functions

This code creates the runqpbox4 helper function.

```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```

This code creates the `runqpbox4prec` helper function.

```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```