Большая разреженная квадратичная программа, основанная на проблеме

Этот пример показывает значение использования разреженной арифметики, когда у вас есть разреженная проблема. Матрица имеет строки n, где вы выбираете n, чтобы быть большим значением и несколькими ненулевыми диагональными полосами. Полная матрица размера, n-by-n может израсходовать всю доступную память, но разреженная матрица не представляет проблемы.

Проблема состоит в том, чтобы минимизировать x'*H*x/2 + f'*x, подвергающийся

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

где f = [-1;-2;-3;...;-n]. H является разреженной симметричной ленточной матрицей.

Создайте разреженную квадратичную матрицу

Создайте симметричную циркулянтную матрицу H на основе сдвигов векторного [3,6,2,14,2,6,3], с 14 находиться на основной диагонали. Имейте матрицу быть n-by-n, где 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);

Просмотрите структуру разреженной матрицы.

spy(H)

Создайте переменные оптимизации и проблему

Создайте переменную x оптимизации и проблему qprob.

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

Создайте целевую функцию и ограничения. Поместите цель и ограничения в qprob.

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

Решите проблему

Решите проблему квадратичного программирования с помощью 'алгоритма interior-point-convex' по умолчанию и линейной алгебры для разреженных матриц. Чтобы помешать решателю останавливаться преждевременно, установите опцию StepTolerance на 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>

Исследуйте решение

Просмотрите значение целевой функции, количество итераций и множитель Лагранжа, сопоставленный с линейным ограничением неравенства.

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.

Оцените ограничение, чтобы видеть, что решение находится на контуре.

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.

Сумма компонентов решения является нулем к в допусках.

Решение x имеет три области: начальный фрагмент, итоговый фрагмент и приблизительно линейный фрагмент по большей части решения. Постройте эти три области.

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

Документация