Main Content

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 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 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 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 value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance. <stopping criteria details>

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 7. 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 7.599738e-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)')