Quadratic Programming Examples

On this page…
Example: Quadratic Minimization with Bound Constraints Example: Quadratic Minimization with a Dense but Structured Hessian Example: Large Sparse Quadratic Program with Interior Point Algorithm

Example: Quadratic Minimization with Bound Constraints

To minimize a large-scale quadratic with upper and lower bounds, you can use the quadprog function.

The problem stored in the MAT-file qpbox1.mat is a positive definite quadratic, and the Hessian matrix H is tridiagonal, subject to upper (ub) and lower (lb) bounds.

Step 1: Load the Hessian and define f, lb, and ub.

load qpbox1   % Get H
lb = zeros(400,1); lb(400) = -inf;
ub = 0.9*ones(400,1); ub(400) = inf;
f = zeros(400,1); f([1 400]) = -2;

Step 2: Call a quadratic minimization routine with a starting point xstart.

xstart = 0.5*ones(400,1);
[x,fval,exitflag,output] = ... 
        quadprog(H,f,[],[],[],[],lb,ub,xstart);

Looking at the resulting values of exitflag and output,

exitflag,output

exitflag =
     3

output = 
          algorithm: 'trust-region-reflective'
         iterations: 19
    constrviolation: 0
      firstorderopt: 1.3674e-005
       cgiterations: 1644
            message: [1x206 char]

You can see that while convergence occurred in 20 iterations, the high number of CG iterations indicates that the cost of the linear system solve is high. In light of this cost, one strategy would be to limit the number of CG iterations per optimization iteration. The default number is the dimension of the problem divided by two, 200 for this problem. Suppose you limit it to 50 using the MaxPCGIter flag in options:

options = optimset('MaxPCGIter',50);
[x,fval,exitflag,output] = ... 
        quadprog(H,f,[],[],[],[],lb,ub,xstart,options);

This time convergence still occurs and the total number of CG iterations (1547) has dropped:

exitflag,output

exitflag =
     3

output = 
          algorithm: 'trust-region-reflective'
         iterations: 36
    constrviolation: 0
      firstorderopt: 2.3821e-005
       cgiterations: 1547
            message: [1x206 char]

A second strategy would be to use a direct solver at each iteration by setting the PrecondBandWidth option to inf:

options = optimset('PrecondBandWidth',inf);
[x,fval,exitflag,output] = ... 
        quadprog(H,f,[],[],[],[],lb,ub,xstart,options);

Now the number of iterations has dropped to 10:

exitflag,output

exitflag =
     3

output = 
          algorithm: 'trust-region-reflective'
         iterations: 10
    constrviolation: 0
      firstorderopt: 6.0684e-009
       cgiterations: 0
            message: [1x206 char]

Using a direct solver at each iteration usually causes the number of iterations to decrease, but often takes more time per iteration. For this problem, the tradeoff is beneficial, as the time for quadprog to solve the problem decreases by a factor of 10.

Back to Top

Example: Quadratic Minimization with a Dense but Structured Hessian

The quadprog large-scale 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 medium-scale problems and for large-scale 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.

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

quadprog(B,f,[],[],[],[],l,u,xstart,options)

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 form

  W = qpbox4mult(Hinfo,Y,...)

  The first two arguments Hinfo and Y are required.

• Loads the problem parameters from qpbox4.mat.

• Uses optimset 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.

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  
% subfunction

% Choose the HessMult option
options = optimset('HessMult',mtxmpy);

% Pass B to qpbox4mult via the Hinfo 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 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: [1x206 char]

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

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 
subfunction

% Choose the HessMult option
% Override the TolPCG option
options = optimset('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 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

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

Note   Decreasing TolPcg too much can substantially increase the number of PCG iterations.

Back to Top

Example: Large Sparse Quadratic Program with Interior Point Algorithm

This example shows the value of using sparse arithmetic when you have a sparse problem. The matrix has n rows, where you choose n. For some large n, the active-set algorithm runs out of memory, but the interior-point-convex algorithm works fine.

The problem is to minimize x'*H*x/2 + f'*x subject to

x(1) + x(2) + ... + x(n) = 0,

where f = [-1;-2;-3;...;-n].

1. Create the parameter n and the utility matrix T. The matrix T is a sparse circulant matrix that is simply a helper for creating the sparse positive-definite quadratic matrix H.

   n = 30000; % Adjust n to a large value
   T = spalloc(n,n,n); % make a sparse circulant matrix
   r = 1:n-1;
   for m = r
       T(m,m+1)=1;
   end
   T(n,1) = 1;

2. Create a sparse vector v. Then create the matrix H by shifted versions of v*v'. The matrix T creates shifts of v.

   v(n) = 0; v(1) = 1; v(2) = 2; v(4) = 3;
   v = (sparse(v))';
   % Make a banded type of matrix
   H = spalloc(n,n,7*n);
   r = 1:n;
   for m = r
       H = H + v*v';
       v = T*v;
   end

3. Take a look at the structure of H:

   spy(H)

   

4. Create the problem vector f and linear constraint.

   f = -r; % linear term
   A = ones(1,n); b = 0;

5. Solve the quadratic programming problem with the interior-point-convex algorithm.

   options = optimset('Algorithm','interior-point-convex');
   [x,fval,exitflag,output,lambda] = ...
       quadprog(H,f,A,b,[],[],[],[],[],options);
   
   Minimum found that satisfies the constraints.
   Optimization completed because the objective function is 
   non-decreasing in feasible directions, to within the selected
   value of the function tolerance, and constraints are satisfied 
   to within the selected value of the constraint tolerance.

6. View the solution value, output structure, and Lagrange multiplier:

   fval,output,lambda
   
   fval =
    -3.1331e+010
   
   output = 
               message: [1x912 char]
             algorithm: 'interior-point-convex'
         firstorderopt: 1.3691e-004
       constrviolation: 9.1268e-009
            iterations: 5
          cgiterations: []
   
   lambda = 
       ineqlin: 1.5000e+004
         eqlin: [0x1 double]
         lower: [30000x1 double]
         upper: [30000x1 double]

   Since there are no lower bounds or upper bounds, all the values in lambda.lower and lambda.upper are 0. The inequality constraint is active, since lambda.ineqlin is nonzero.

7. Notice that quadprog with the active-set algorithm fails with an out-of-memory error:

   options = optimset('Algorithm','active-set');
   [x fval] = quadprog(H,f,A,b,[],[],[],[],[],options);
   
   Warning: Cannot use sparse matrices with active-set algorithm: 
       converting to full. 
   > In quadprog at 377
   Error using eye
   Out of memory. Type HELP MEMORY for your options.
   
   Error in qpsub at 224
           Q = eye(numberOfVariables,numberOfVariables);
   
   Error in quadprog at 429
      [X,lambdaqp,exitflag,output,~,~,msg]= ...

Back to Top

Quadratic Programming Algorithms

Binary Integer Programming Algorithms

© 1984-2012- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS