%% 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);

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

% 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 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)')

% Copyright 2018 The MathWorks, Inc.