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

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

Создайте линейное ограничение и цель

Линейное ограничение состоит в том, что сумма элементов решения неположительна. Целевая функция содержит линейный член, выраженный в векторном f.

A = ones(1,n);
b = 0;
f = 1:n;
f = -f;

Решите задачу

Решите задачу квадратичного программирования с помощью 'interior-point-convex' алгоритм. Чтобы помешать решателю останавливаться преждевременно, установите StepTolerance опция к 0.

options = optimoptions(@quadprog,'Algorithm','interior-point-convex','StepTolerance',0);
[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 value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

На многих компьютерах вы не можете создать полный n- n матрица, когда n = 30,000. Таким образом, можно запустить эту проблему только при помощи разреженных матриц.

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

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

fprintf('The objective function value is %d.\nThe number of iterations is %d.\nThe Lagrange multiplier is %d.\n',...
    fval,output.iterations,lambda.ineqlin)

The objective function value is -3.133073e+10.
The number of iterations is 7.
The Lagrange multiplier is 1.500050e+04.

Поскольку нет никаких нижних границ, верхних границ или линейных ограничений равенства, единственным значимым множителем Лагранжа является lambda.ineqlin. Поскольку lambda.ineqlin является ненулевым, можно сказать, что ограничение неравенства активно. Оцените ограничение, чтобы видеть, что решение находится на контуре.

fprintf('The linear inequality constraint A*x has value %d.\n',A*x)

The linear inequality constraint A*x has value 9.150244e-08.

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

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

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