This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

Large Sparse Quadratic Program, Problem-Based

This example shows the value of using sparse arithmetic when you have a sparse problem. The matrix has n rows, where you choose n to be a large value, and a few nonzero diagonal bands. A full matrix of size n-by-n can use up all available memory, but a sparse matrix presents no problem.

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]. H is a sparse symmetric banded matrix.

Create Sparse Quadratic Matrix

Create a symmetric circulant matrix H based on shifts of the vector [3,6,2,14,2,6,3], with 14 being on the main diagonal. Have the matrix be n-by-n, where n = 30,000.

n = 3e4;
H2 = speye(n);
H = 3*circshift(H2,-3,2) + 6*circshift(H2,-2,2) + 2*circshift(H2,-1,2)...
    + 14*H2 + 2*circshift(H2,1,2) + 6*circshift(H2,2,2) + 3*circshift(H2,3,2);

View the sparse matrix structure.

spy(H)

Create Optimization Variables and Problem

Create an optimization variable x and problem qprob.

x = optimvar('x',n);
qprob = optimproblem;

Create the objective function and constraints. Place the objective and constraints into qprob.

f = 1:n;
obj = 1/2*x'*H*x - f*x;
qprob.Objective = obj;
cons = sum(x) <= 0;
qprob.Constraints = cons;

Solve Problem

Solve the quadratic programming problem using the default 'interior-point-convex' algorithm and sparse linear algebra. To keep the solver from stopping prematurely, set the StepTolerance option to 0.

options = optimoptions('quadprog','Algorithm','interior-point-convex',...
    'LinearSolver','sparse','StepTolerance',0);
[sol,fval,exitflag,output,lambda] = solve(qprob,'Options',options);
Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

Examine Solution

View the objective function value, number of iterations, and Lagrange multiplier associated with the linear inequality constraint.

fprintf('The objective function value is %d.\nThe number of iterations is %d.\nThe Lagrange multiplier is %d.\n',...
    fval,output.iterations,lambda.Constraints)
The objective function value is -3.133073e+10.
The number of iterations is 6.
The Lagrange multiplier is 1.500050e+04.

Evaluate the constraint to see that the solution is on the boundary.

fprintf('The linear inequality constraint sum(x) has value %d.\n',sum(sol.x))
The linear inequality constraint sum(x) has value 9.726136e-09.

The sum of the solution components is zero to within tolerances.

The solution x has three regions: an initial portion, a final portion, and an approximately linear portion over most of the solution. Plot the three regions.

subplot(3,1,1)
plot(sol.x(1:60))
title('x(1) Through x(60)')
subplot(3,1,2)
plot(sol.x(61:n-60))
title('x(61) Through x(n-60)')
subplot(3,1,3)
plot(sol.x(n-59:n))
title('x(n-59) Through x(n)')

Related Topics