Квадратичные алгоритмы программирования

Определение квадратичного программирования

Квадратичное программирование - это задача нахождения вектора x, которая минимизирует квадратичную функцию, возможно, удовлетворяющую линейным ограничениям:

\min_{x} \frac{1}{2} x^{T} H x + c^{T} x

(1)

таким образом, что <reservedrangesplaceholder6> ≤ <reservedrangesplaceholder5> , Aeq·x = beq, <reservedrangesplaceholder2> ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>.

`interior-point-convex` `quadprog` Алгоритм

The interior-point-convex алгоритм выполняет следующие шаги:

Предварительное решение/Постсольве
Сгенерируйте начальную точку
Предиктор-корректор
Условия остановки
Обнаружение недопустимости

Примечание

Алгоритм имеет два пути кода. Это принимает один, когда H матрицы Гессия является обыкновенной (полной) матрицей двойников, и он принимает другой, когда H является разреженной матрицей. Для получения дополнительной информации о разреженном типе данных см. «Разреженные матрицы». Обычно алгоритм быстрее для больших задач, которые имеют относительно мало ненулевых членов, когда вы задаете H как sparse. Точно так же алгоритм быстрее для малых или относительно плотных задач, когда вы задаете H как full.

Предварительное решение/Постсольве

Алгоритм сначала пытается упростить задачу, удаляя избыточности и упрощая ограничения. Задачи, выполненные во время предварительного решения, могут включать следующее:

Проверьте, имеют ли какие-либо переменные равные верхнюю и нижнюю границы. Если это так, проверьте допустимость, а затем исправьте и удалите переменные.
Проверяйте, включает ли любое линейное ограничение неравенства только одну переменную. Если это так, проверьте допустимость, а затем измените линейное ограничение на ограничение.
Проверяйте, включает ли любое линейное ограничение равенства только одну переменную. Если это так, проверьте допустимость, а затем исправьте и удалите переменную.
Проверьте, имеет ли какая-либо линейная матрица ограничений нулевые строки. Если это так, проверьте допустимость, а затем удалите строки.
Определите, являются ли границы и линейные ограничения допустимыми.
Проверяйте, появляются ли какие-либо переменные только как линейные условия в целевой функции и не появляются ли они в каких-либо линейных ограничениях. Если это так, проверяйте допустимость и граничность, а затем фиксируйте переменные в соответствующих границах.
Измените любые линейные ограничения неравенства на линейные ограничения равенства путем добавления переменных slack.

Если алгоритм обнаруживает недопустимую или неограниченную проблему, он останавливается и выдает соответствующее выходное сообщение.

Алгоритм может прийти в одну допустимую точку, которая представляет решение.

Если алгоритм не обнаруживает недопустимую или неограниченную задачу в шаге prevolve, и если предварительное решение не выдало решение, алгоритм переходит к своим следующим шагам. После достижения критерия остановки алгоритм восстанавливает исходную задачу, отменяя любые предварительные преобразования. Этот заключительный шаг является шагом после решения.

Для получения дополнительной информации смотрите Gould and Toint [63].

Сгенерируйте начальную точку

Начальная точка x0 для алгоритма является:

Инициализация x0 на ones(n,1), где n - количество строк в H.
Для компонентов, которые имеют обе верхние границы ub и нижнюю границу lb, если компонент x0 не находится строго внутри границ, для компонента задано значение (ub + lb)/2.
Для компонентов, которые имеют только одну привязку, измените компонент, если необходимо, чтобы он находился строго внутри привязки.
Сделайте шаг предиктора (см. Predictor-Corrector), с незначительными коррекциями для выполнимости, а не полный шаг предиктора-корректора. Это помещает начальную точку ближе к central path, не влекя за собой накладные расходы на полный шаг предиктор-корректор. Для получения дополнительной информации о центральном пути смотрите Nocedal и Wright [7], страница 397.

Предиктор-корректор

Разреженные и полные промежуточна точка выпуклые алгоритмы различаются в основном фазой предиктор-корректор. Алгоритмы похожи, но отличаются некоторыми деталями. Описание базового алгоритма смотрите в Mehrotra [47].

Алгоритмы начинаются с превращения линейных неравенств Ax < = b в неравенства вида Ax > = b умножением A и b на -1. Это не имеет никакого значения для решения, но делает проблему той же формы, что и в некоторой литературе.

Разреженный предиктор-корректор
Полный предиктор-корректор

Разреженный предиктор-корректор. Подобно fmincon алгоритм внутренней точки, разреженный interior-point-convex алгоритм пытается найти точку, в которой держатся условия Каруша-Куна-Такера (KKT). Для квадратичной задачи программирования, описанной в определении квадратичного программирования, эти условия следующие:

$\begin{matrix} H x + c - A_{e q}^{T} y - {\bar{A}}^{T} z = 0 \\ \bar{A} x - \bar{b} - s = 0 \\ A_{e q} x - b_{e q} = 0 \\ s_{i} z_{i} = 0, i = 1, 2, ..., m \\ s \geq 0 \\ z \geq 0. \end{matrix}$

Здесь

$\bar{A}$ - расширенная линейная матрица неравенства, которая включает ограничения, записанные как линейные неравенства. $\bar{b}$ - соответствующий вектор линейного неравенства, включая ограничения.
s - вектор ослаблений, которые преобразуют ограничения неравенства в равенства. s имеет m длины, количество линейных неравенств и границ.
z - вектор множителей Лагранжа, соответствующих s.
y - вектор множителей Лагранжа, сопоставленный с ограничениями равенства.

Алгоритм сначала предсказывает шаг из формулы Ньютона-Рафсона, затем вычисляет шаг корректора. Корректор пытается лучше применить нелинейное ограничение _sizi = 0.

Определения для шага предиктора:

_rd, двойная невязка:
$r_{d} = H x + c - A_{e q}^{T} y - {\bar{A}}^{T} z .$
_req основное ограничение равенства невязки:
$r_{e q} = A_{e q} x - b_{e q} .$
_rineq, основная невязка ограничения неравенства, который включает ограничения и ослабления:
$r_{i n e q} = \bar{A} x - \bar{b} - s .$
_rsz комплементарность невязки:
_rsz = Sz.
S - диагональная матрица членов slack, z - матрица столбцов множителей Лагранжа.
_rc средняя комплементарность:
$r_{c} = \frac{s^{T} z}{m} .$

В шаге Ньютона изменения в x, s, y и z заданы:

(\begin{matrix} H & 0 & - A_{e q}^{T} & - {\bar{A}}^{T} \\ A_{e q} & 0 & 0 & 0 \\ \bar{A} & - I & 0 & 0 \\ 0 & Z & 0 & S \end{matrix}) (\begin{matrix} Δ x \\ Δ s \\ Δ y \\ Δ z \end{matrix}) = - (\begin{matrix} r_{d} \\ r_{e q} \\ r_{i n e q} \\ r_{s z} \end{matrix}) .

(2)

Однако полный шаг Ньютона может быть недопустимым, из-за ограничений позитивности на s и z. Поэтому, quadprog сокращает шаг, при необходимости, для поддержания позитива.

Кроме того, чтобы сохранить «центрированное» положение во внутреннем пространстве, вместо того, чтобы пытаться решить _sizi = 0, алгоритм принимает положительное σ параметра и пытается решить

_sizi = _σrc.

quadprog заменяет _rsz в уравнении шага Ньютона с _rsz + Δ <reservedrangesplaceholder2> Δ <reservedrangesplaceholder1> - <reservedrangesplaceholder0> 1, где 1 вектор таковых. Также,quadprog переупорядочивает уравнения Ньютона, чтобы получить симметричную, более численно стабильную систему для вычисления шага предиктора.

После вычисления исправленного шага Ньютона алгоритм выполняет больше вычислений, чтобы получить как более длинный текущий шаг, так и подготовиться к лучшим последующим шагам. Эти многочисленные вычисления коррекции могут улучшить как эффективность, так и робастность. Для получения дополнительной информации смотрите Gondzio [4].

Полный предиктор-корректор. Полный алгоритм предиктор-корректор не объединяет границы в линейные ограничения, поэтому у него есть другой набор переменных slack, соответствующих границам. Алгоритм смещает нижние границы в нуль. И, если существует только одна граница для переменной, алгоритм превращает ее в нижнюю границу в нуле, отрицая неравенство верхней границы.

$\bar{A}$ - расширенная линейная матрица, которая включает как линейные неравенства, так и линейные равенства. $\bar{b}$ - соответствующий линейный вектор равенства. $\bar{A}$ также включает условия расширения вектора x с переменными slack s которые поворачивают ограничения неравенства к ограничениям равенства:

$\bar{A} x = (\begin{matrix} A_{e q} & 0 \\ A & I \end{matrix}) (\begin{matrix} x_{0} \\ s \end{matrix}),$

где x 0 означает исходный вектор x.

Условия KKT:

\begin{matrix} H x + c - {\bar{A}}^{T} y - v + w = 0 \\ \bar{A} x = \bar{b} \\ x + t = u \\ v_{i} x_{i} = 0, i = 1, 2, ..., m \\ w_{i} t_{i} = 0, i = 1, 2, ..., n \\ x, v, w, t \geq 0. \end{matrix}

(3)

Чтобы найти x решения, переменные slack и двойственные переменные к Уравнению 3, алгоритм в основном рассматривает шаг Ньютона-Рафсона:

(\begin{matrix} H & - {\bar{A}}^{T} & 0 & - I & I \\ \bar{A} & 0 & 0 & 0 & 0 \\ - I & 0 & - I & 0 & 0 \\ V & 0 & 0 & X & 0 \\ 0 & 0 & W & 0 & T \end{matrix}) (\begin{matrix} Δ x \\ Δ y \\ Δ t \\ \begin{array}{l} Δ v \\ Δ w \end{array} \end{matrix}) = - (\begin{matrix} H x + c - {\bar{A}}^{T} y - v + w \\ \bar{A} x - \bar{b} \\ u - x - t \\ \begin{array}{l} V X \\ W T \end{array} \end{matrix}) = - (\begin{matrix} r_{d} \\ r_{p} \\ r_{u b} \\ \begin{array}{l} r_{v x} \\ r_{w t} \end{array} \end{matrix}),

(4)

где X, V, W, и T являются диагональными матрицами, соответствующими векторам x, v, w, и t соответственно. Векторы невязок на крайней правой стороне уравнения:

_rd, двойной остаточный
_rp, основная невязка
_rub верхняя граничная невязка
_rvx нижняя граница комплементарности невязки
_rwt верхняя граница комплементарности невязки

Алгоритм решает Уравнение 4, сначала преобразовав его в симметричную матричную форму

(\begin{matrix} - D & {\bar{A}}^{T} \\ \bar{A} & 0 \end{matrix}) (\begin{matrix} Δ x \\ Δ y \end{matrix}) = - (\begin{matrix} R \\ r_{p} \end{matrix}),

(5)

где

$\begin{array}{l} D = H + X^{- 1} V + T^{- 1} W \\ R = - r_{d} - X^{- 1} r_{v x} + T^{- 1} r_{w t} + T^{- 1} W r_{u b} . \end{array}$

Все обратные матрицы в определениях D и R просты в вычислении, потому что матрицы диагональны.

Чтобы вывести Уравнение 5 из Уравнения 4, заметьте, что вторая строка Уравнения 5 такая же, как и вторая строка матрицы Уравнения 4. Первая строка Уравнения 5 происходит от решения двух последних строк Уравнения 4 для В v и В w, а затем решения для В t.

Чтобы решить Уравнение 5, алгоритм следует существенным элементам Altman и Gondzio [1]. Алгоритм решает симметричную систему разложением LDL. Как указывали такие авторы, как Вандербей и Карпентер [2], это разложение численно стабильно без какого-либо поворота, поэтому может быть быстрым.

Полное quadprog алгоритм предиктор-корректор в значительной степени тот же, что и в linprog 'interior-point' алгоритм, но включает и квадратичные условия. См. «Предиктор-Корректор».

Ссылки

[1] Альтман, Анна и Дж. Гондзио. Регуляризованные симметричные неопределенные системы во внутренних точечных методах для линейной и квадратичной оптимизации. Методы и программное обеспечение оптимизации, 1999. Доступно для загрузки здесь.

[2] Вандербей, Р. Дж. и Т. Дж. Карпентер. Симметричные неопределенные системы для внутренних методов точки. Математическое программирование 58, 1993. стр 1–32. Доступно для загрузки здесь.

Условия остановки

Алгоритм предиктора-корректора итерации пока не достигает точки, которая является допустимой (удовлетворяет ограничениям в пределах допусков) и где относительные размеры шага малы. В частности, задайте

$ρ = \max (1, ‖ H ‖, ‖ \bar{A} ‖, ‖ A_{e q} ‖, ‖ c ‖, ‖ \bar{b} ‖, ‖ b_{e q} ‖)$

Алгоритм останавливается, когда все эти условия удовлетворены:

$\begin{matrix} {‖ r_{p} ‖}_{1} + {‖ r_{u b} ‖}_{1} \leq ρ TolCon \\ {‖ r_{d} ‖}_{\infty} \leq ρ TolFun \\ r_{c} \leq TolFun, \end{matrix}$

где

$r_{c} = \max_{i} (\min (| x_{i} v_{i} |, | x_{i} |, | v_{i} |), \min (| t_{i} w_{i} |, | t_{i} |, | w_{i} |)) .$

_rc по существу измеряет размер невязок комплементарности xv и tw, которые являются каждыми нулевыми векторами в решении.

Обнаружение недопустимости

quadprog вычисляет функцию заслуг, φ при каждой итерации. Функция заслуг является мерой выполнимости. quadprog останавливается, если функция заслуг становится слишком большой. В этом случае, quadprog объявляет проблему недопустимой.

Функция заслуг связана с условиями KKT для задачи - см. Предсказатель-Корректор. Используйте следующие определения:

$\begin{matrix} ρ = \max (1, ‖ H ‖, ‖ \bar{A} ‖, ‖ A_{e q} ‖, ‖ c ‖, ‖ \bar{b} ‖, ‖ b_{e q} ‖) \\ r_{eq} = A_{eq} x - b_{eq} \\ r_{ineq} = \bar{A} x - \bar{b} - s \\ r_{d} = H x + c + A_{eq}^{T} λ_{eq} + {\bar{A}}^{T} {\bar{λ}}_{ineq} \\ g = \frac{1}{2} x^{T} H x + c^{T} x - {\bar{b}}^{T} {\bar{λ}}_{ineq} - b_{eq}^{T} λ_{eq} . \end{matrix}$

Обозначение $\bar{A}$ и $\bar{b}$ означает линейные коэффициенты неравенства, дополненные терминами, чтобы представлять границы для разреженного алгоритма. Обозначение ${\bar{λ}}_{ineq}$ аналогично представляет множители Лагранжа для линейных ограничений неравенства, включая связанные ограничения. Это называлось z в Predictor-Corrector, и $λ_{eq}$ назывался y.

Функция заслуг φ

$\frac{1}{ρ} (\max ({‖ r_{eq} ‖}_{\infty}, {‖ r_{ineq} ‖}_{\infty}, {‖ r_{d} ‖}_{\infty}) + g) .$

Если эта функция заслуг становится слишком большой, quadprog объявляет задачу недопустимой и останавливается с выходным флагом -2.

`trust-region-reflective` `quadprog` Алгоритм

Многие методы, используемые в решателях Optimization Toolbox™, основаны на доверительных областях, простой, но мощной концепции в оптимизации.

Чтобы понять подход доверительной области к оптимизации, рассмотрим задачу минимизации без ограничений, минимизируйте f (x), где функция принимает векторные аргументы и возвращает скаляры. Предположим, что вы находитесь в точке, x в n -пространство, и вы хотите улучшить, т.е. переместиться в точку с более низким значением функции. Основная идея состоит в том, чтобы аппроксимировать f с более простой функциональной q, которая разумно отражает поведение функциональной f в окрестности N вокруг точечной x. Этот район - доверительная область. Пробный шаг s вычисляется путем минимизации (или приблизительно минимизации) по N. Это подпрограмма доверительной области,

\min_{s} {q (s), s \in N} .

(6)

Текущая точка обновляется следующим образом x + s, если f ( x + s ) < f (x); в противном случае текущая точка остается неизменной и N, область доверия, сокращается, и расчет пробного шага повторяется.

Ключевыми вопросами в определении конкретного подхода доверительной области к минимизации f (x) являются то, как выбрать и вычислить q приближения (заданную в текущей точке x), как выбрать и изменить N доверительной области и как точно решить подблем доверительной области. Этот раздел посвящен проблеме без ограничений. Более поздние разделы обсуждают дополнительные осложнения из-за наличия ограничений на переменные.

В стандартном методе доверительной области ([48]) квадратичное q приближения определяется первыми двумя слагаемыми приближения Тейлора к F при x; окрестность N обычно сферическая или эллипсоидальная по форме. Математически подпрограмма доверительной области обычно указывается

\min {\frac{1}{2} s^{T} H s + s^{T} g таким , что ‖ D s ‖ \leq Δ},

(7)

где g - градиент f в текущей точке x, H - матрица Мешковины (симметричная матрица вторых производных), D - диагональная матрица масштабирования, Δ - положительная скалярная величина и ∥. ∥ - 2-норма. Хорошие алгоритмы существуют для решения уравнения 7 (см. [48]); такие алгоритмы обычно включают в себя расчет всех собственных значений H и процесс Ньютона, примененный к светскому уравнению

$\frac{1}{Δ} - \frac{1}{‖ s ‖} = 0.$

Такие алгоритмы обеспечивают точное решение уравнения 7. Однако они требуют времени, пропорционального нескольким факторизациям H. Поэтому для масштабных проблем нужен другой подход. В литературе было предложено несколько аппроксимационных и эвристических стратегий, основанных на уравнении 7 ([42] и [50]). Подход приближения, используемый в решателях Optimization Toolbox, состоит в том, чтобы ограничить подблиму доверительной области двумерной подпространственной S ([39] и [42]). После вычисления S подпространства работа по решению Уравнения 7 тривиальна, даже если необходима полная собственная/собственная информация (поскольку в подпространстве задача является только двумерной). Доминирующая работа теперь сместилась на определение подпространства.

Двумерный подпространственный S определяется с помощью предварительно обусловленного сопряженного градиентного процесса, описанного ниже. Решатель определяет S как линейное пространство, охватываемое s 1 и _s 2, где s 1 находится в направлении градиента g, и s 2 является либо приблизительным направлением Ньютона, т.е. решением

H \cdot s_{2} = - g,

(8)

или направление отрицательной кривизны,

s_{2}^{T} \cdot H \cdot s_{2} < 0.

(9)

Философия этого выбора S состоит в том, чтобы форсировать глобальное сходимость (через направление наискорейшего спуска или направление отрицательной кривизны) и достичь быстрого локального сходимости (через шаг Ньютона, когда он существует).

Эскиз без ограничений минимизации с использованием идей доверительной области теперь легко дать:

Сформулируйте двумерную подпрограмму доверительной области.
Решите уравнение 7, чтобы определить пробный шаг s.
Если f (x + s) <f (<reservedrangesplaceholder3>), то x = x + s.
Отрегулируйте

Эти четыре шага повторяются до сходимости. Размерность Β доверительной области регулируется согласно стандартным правилам. В частности, это уменьшено, если шаг испытания не принят, т.е. f (x + s) ≥ <reservedrangesplaceholder1> (<reservedrangesplaceholder0>). Для обсуждения этого аспекта см. [46] и [49].

Решатели Optimization Toolbox обрабатывают несколько важных особых случаев f со специализированными функциями: нелинейные методы наименьших квадратов, квадратичные функции и линейные методы наименьших квадратов. Однако базовые алгоритмические идеи те же, что и в общем случае. Эти особые случаи рассматриваются в последующих разделах.

Метод доверительной области подпространства используется для определения направления поиска. Однако вместо ограничения шага (возможно) одним шагом отражения, как в случае нелинейной минимизации, при каждой итерации проводится кусочно-отражающий поиск линии. Для получения дополнительной информации о поиске по линии см. раздел [45].

Предварительно обусловленный сопряженный градиентный метод

Популярным способом решения больших, симметричных, положительно определенных систем линейных уравнений Hp = - g является метод Предварительно обусловленных сопряженных градиентов (PCG). Этот итеративный подход требует способности вычислять матрично-векторные продукты вида H·v где v является произвольным вектором. Симметричная положительная определенная матричная M является предварительным условием для H. То есть M = C², где C^–1HC^–1 является хорошо обусловленной матрицей или матрицей с кластеризованными собственными значениями.

В контексте минимизации можно предположить, что H матрицы Гессия симметрична. Однако H гарантированно будет положительно определено только в окрестности сильного минимизатора. Алгоритм PCG выходит, когда он сталкивается с направлением отрицательной (или нулевой) кривизны, то есть d^THd ≤ 0. Выходное направление PCG p является либо направлением отрицательной кривизны, либо приблизительным решением системы Ньютона Hp = - g. В любом случае p помогает задать двумерный подпространство, используемый в подходе доверительной области, обсуждаемом в Методах доверительной области для нелинейной минимизации.

Линейные ограничения равенства

Линейные ограничения усложняют ситуацию, описанную для минимизации без ограничений. Однако основные идеи, описанные выше, могут осуществляться чистым и эффективным способом. Методы доверительной области в решателях Optimization Toolbox генерируют строго допустимые итерации.

Общая задача минимизации с ограничениями линейного равенства может быть записана

\min {f (x) such that A x = b},

(10)

где A m n матрицей (m ≤ <reservedrangesplaceholder1>). Некоторые решатели Optimization Toolbox предварительно обрабатывают A, чтобы удалить строгие линейные зависимости, используя метод, основанный на LU-факторизации A^T [46]. Здесь A принято рангом m.

Метод, используемый для решения уравнения 10, отличается от подхода без ограничений двумя значительными способами. Во-первых, начальная допустимая точка x 0 вычисляется, используя разреженный шаг наименьших квадратов, так что Ax 0 = b. Во-вторых, Алгоритм PCG заменяется на Сокращенные Предварительно Обусловленные Сопряженные Градиенты (RPCG), см. [46], в порядок чтобы вычислить приблизительный уменьшенный шаг Ньютона (или направление отрицательной кривизны в нулевом пространстве A). Ключевой шаг линейной алгебры включает в себя решение систем вида

[\begin{matrix} C & {\tilde{A}}^{T} \\ \tilde{A} & 0 \end{matrix}] [\begin{matrix} s \\ t \end{matrix}] = [\begin{matrix} r \\ 0 \end{matrix}],

(11)

где $\tilde{A}$ аппроксимирует A (маленькие ненулевые A устанавливаются в нули при условии, что ранг не потерян) и C является разреженным симметричным положительно-определенным приближением к H, т.е., C = H. Для получения дополнительной информации см. раздел [46].

Прямоугольные ограничения

Задача с ограничениями в прямоугольнике имеет вид

\min {f (x) such that l \leq x \leq u},

(12)

где l - вектор нижних границ, а u - вектор верхних границ. Некоторые (или все) компоненты l могут быть равны - ∞ и некоторые (или все) компоненты u могут быть равны ∞. Метод генерирует последовательность строго допустимых точек. Два метода используются для поддержания выполнимости при достижении устойчивого поведения сходимости. Во-первых, масштабированный измененный шаг Ньютона заменяет необязательный шаг Ньютона (для определения двумерного подпространства S). Во-вторых, отражения используются для увеличения размера шага.

Масштабированный модифицированный шаг Ньютона возникает из-за исследования необходимых условий Куна-Такера для уравнения 12,

{(D (x))}^{- 2} g = 0,

(13)

где

$D (x) = diag ({| v_{k} |}^{- 1 / 2}),$

и вектор v (<reservedrangesplaceholder2>) определен ниже для каждого 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>:

Если _gi < 0 и _ui < ∞ то _vi = _xi - _ui
Если _gi ≥ 0 и _li > - ∞ то _vi = _xi - _li
Если _gi < 0 и _ui = ∞ то _vi = -1
Если _gi ≥ 0 и _li = - ∞ то _vi = 1

Нелинейная система Уравнение 13 не везде дифференцируема. Недифференцируемость возникает, когда _vi = 0. Можно избегать таких точек, сохраняя строгую допустимость, т.е. ограничивая l < x < u.

Масштабированный модифицированный шаг Ньютона _sk для нелинейной системы уравнений, заданных уравнением 13, задан как решение линейной системы

\hat{M} D s^{N} = - \hat{g}

(14)

на k-й итерации, где

\hat{g} = D^{- 1} g = diag ({| v |}^{1 / 2}) g,

(15)

\hat{M} = D^{- 1} H D^{- 1} + diag (g) J^{v} .

(16)

Здесь J^v играет роль якобиана |<reservedrangesplaceholder0>|. Каждая диагональ компонент матрицы диагонали J^v равен 0, -1 или 1. Если все компоненты l и u конечны, J^v = diag (знак (g)). В точке, где gi = 0, vi может не быть дифференцируемым. $J_{i i}^{v} = 0$ определяется в такой точке. Недифференцируемость этого типа не является причиной для беспокойства, потому что для такого компонента не важно, какое значение принимает _vi. Кроме того, _{|<reservedrangesplaceholder2>|} все еще будет прерывистым в этой _точке, но |<reservedrangesplaceholder1>|·<reservedrangesplaceholder0> функции непрерывна.

Во-вторых, отражения используются для увеличения размера шага. Шаг (одинарный) отражения определяется следующим образом. Учитывая p шага, которая пересекает ограничение привязки, рассмотрим первое ограничение привязки, пересекаемое p; предположим, что это i-е ограничение (i-е верхнее или i-е нижнее ограничение). Затем шаг отражения p^R = p кроме i-го компонента, где p^R_i = - _pi.

`active-set` `quadprog` Алгоритм

После завершения предварительного решения, active-set алгоритм протекает в две фазы.

Фаза 1 - Получение допустимой точки относительно всех ограничений.
Фаза 2 - итерационно понижает целевую функцию с сохранением списка активных ограничений и сохранением выполнимости в каждой итерации.

The active-set стратегия (также известная как метод проекции) аналогична стратегии Gill et al., описанной в [18] и [17].

Шаг предварительного решения

Проверьте, имеют ли какие-либо переменные равные верхнюю и нижнюю границы. Если это так, проверьте допустимость, а затем исправьте и удалите переменные.
Проверяйте, включает ли любое линейное ограничение неравенства только одну переменную. Если это так, проверьте допустимость, а затем измените линейное ограничение на ограничение.
Проверяйте, включает ли любое линейное ограничение равенства только одну переменную. Если это так, проверьте допустимость, а затем исправьте и удалите переменную.
Проверьте, имеет ли какая-либо линейная матрица ограничений нулевые строки. Если это так, проверьте допустимость, а затем удалите строки.
Определите, являются ли границы и линейные ограничения допустимыми.
Проверяйте, появляются ли какие-либо переменные только как линейные условия в целевой функции и не появляются ли они в каких-либо линейных ограничениях. Если это так, проверяйте допустимость и граничность, а затем фиксируйте переменные в соответствующих границах.
Измените любые линейные ограничения неравенства на линейные ограничения равенства путем добавления переменных slack.

Алгоритм может прийти в одну допустимую точку, которая представляет решение.

Для получения дополнительной информации смотрите Gould and Toint [63].

Алгоритм фазы 1

В Фазе 1 алгоритм пытается найти точку x который удовлетворяет всем ограничениям без учета целевой функции. quadprog запускает алгоритм Фаза 1 только в том случае, если заданная начальная точка x0 недопустимо.

Чтобы начать, алгоритм пытается найти точку, которая допустима относительно всех ограничений равенства, таких как X = Aeq\beq. Если решения нет x к уравнениям Aeq*x = beq, затем алгоритм останавливается. Если существует решение Xследующий шаг состоит в том, чтобы удовлетворить границам и линейным неравенствам. В случае отсутствия ограничений равенства X = x0, начальная точка.

Начиная с Xалгоритм вычисляет M = max(A*X – b, X - ub, lb – X). Если M <= options.ConstraintTolerance, затем точка X является допустимым, и алгоритм Фазы 1 останавливается.

Если M > options.ConstraintToleranceалгоритм представляет неотрицательную переменную slack γ для вспомогательной задачи линейного программирования

$\min_{x, γ} γ$

таким, что

$\begin{matrix} A x - γ \leq b \\ Aeq x = beq \\ фунт - γ \leq x \leq ub + γ \\ γ \geq - ρ . \end{matrix}$

Вот, ρ есть ConstraintTolerance опция умножена на абсолютное значение самого большого элемента в A и Aeq. Если алгоритм находит γ = 0 и x точки, которая удовлетворяет уравнениям и неравенствам, то x является допустимой точкой Фазы 1. Если нет решения вспомогательной задачи линейного программирования x с γ = 0, то задача Фазы 1 недопустима.

Чтобы решить задачу вспомогательного линейного программирования, алгоритм задает γ 0 = M + 1, устанавливает _x 0 = X, и инициализирует активный набор как фиксированные переменные (если таковые имеются) и все ограничения равенства. Алгоритм переформулирует переменные линейного программирования, p как смещение x от текущей точки x 0, а именно x = _x 0 + p. Алгоритм решает задачу линейного программирования теми же итерациями, что и в Фазе 2, чтобы решить квадратичную задачу программирования с соответствующим образом измененным Гессианом.

Алгоритм фазы 2

С точки зрения переменной d задача заключается в том,

\begin{matrix} \min_{d \in ℜ^{n}} q (d) = \frac{1}{2} d^{T} H d + c^{T} d, \\ A_{i} d = b_{i}, i = 1, ..., m_{e} \\ A_{i} d \leq b_{i}, i = m_{e} + 1, ..., m . \end{matrix}

(17)

Здесь _Ai относится к iвторая строка m -by n матрицы A.

Во время фазы 2 активный набор ${\bar{A}}_{k}$ , которая является оценкой активных ограничений (ограничений на границах ограничений) в точке решения.

Алгоритм обновляется ${\bar{A}}_{k}$ на каждом k итерации формируют базис для _dk направления поиска. Ограничения равенства всегда остаются в активном наборе ${\bar{A}}_{k}$ . _dk направления поиска вычисляется и минимизирует целевую функцию, оставаясь на любых активных границах ограничений. Алгоритм формирует допустимое подпространство для _dk из базиса _Zk, столбцы которого ортогональны оценке активного множества ${\bar{A}}_{k}$ (то есть, ${\bar{A}}_{k} Z_{k} = 0$ ). Поэтому направление поиска, которое формируется из линейного суммирования любой комбинации столбцов _Zk, гарантированно остается на контурах активных ограничений.

Алгоритм формирует матрицу, _Zk из последних m - l столбцов QR-разложения матрицы ${\bar{A}}_{k}^{T}$ , где l количество активных ограничений и l < m. То есть _Zk задается как

Z_{k} = Q [:, l + 1 : m],

(18)

где

$Q^{T} {\bar{A}}_{k}^{T} = [\begin{matrix} R \\ 0 \end{matrix}] .$

После нахождения _Zk алгоритм ищет новое _dk направления поиска, которое минимизирует q (d), где _dk находится в нулевом пространстве активных ограничений. То есть _dk является линейной комбинацией столбцов _Zk : ${\hat{d}}_{k} = Z_{k} p$ для некоторых векторных p.

Просмотр квадратичной функции как функции p путем подстановки _dk, дает

q (p) = \frac{1}{2} p^{T} Z_{k}^{T} H Z_{k} p + c^{T} Z_{k} p .

(19)

Дифференцирование этого уравнения относительно p выражений

\nabla q (p) = Z_{k}^{T} H Z_{k} p + Z_{k}^{T} c .

(20)

∇ q (p) упоминается как проективный градиент квадратичной функции, потому что это градиент, проецируемый в подпространстве, заданном Zk. Термин $Z_{k}^{T} H Z_{k}$ называется спроецированным Гессианом. Принимая спроецированную матрицу Гессия $Z_{k}^{T} H Z_{k}$ положительный полуопределенный, минимум функции q (p) в заданном _Zk подпространстве возникает, когда ∇ q (p ) = 0, что является решением системы линейных уравнений

Z_{k}^{T} H Z_{k} p = - Z_{k}^{T} c .

(21)

Алгоритм затем принимает шаг вида

$x_{k + 1} = x_{k} + α d_{k},$

где

$d_{k} = Z_{k} p .$

Из-за квадратичного характера целевой функции на каждой итерации α существовать только два выбора длины шага. Шаг единства вдоль _dk является точным шагом к минимуму функции, ограниченной ядром пространства ${\bar{A}}_{k}$ . Если алгоритм может сделать такой шаг, не нарушая ограничений, то этот шаг является решением квадратичной программы (Уравнение 18). В противном случае шаг вдоль _dk к ближайшему ограничению меньше, чем единица, и алгоритм включает новое ограничение в активном наборе при следующей итерации. Расстояние до границ ограничений в любом направлении _dk задается

$α = \min_{i \in {1, ..., m}} {\frac{- (A_{i} x_{k} - b_{i})}{A_{i} d_{k}}},$

который задан для ограничений, не входящих в активный набор, и где _dk направления направлено к границе ограничений, то есть, $A_{i} d_{k} > 0, i = 1, ..., m$ .

Когда активное множество включает n независимых ограничений, без расположения минимума, алгоритм вычисляет _λk множителей Лагранжа, которые удовлетворяют несингулярному набору линейных уравнений

{\bar{A}}_{k}^{T} λ_{k} = c + H x_{k} .

(22)

Если все элементы _λk положительны, _xk является оптимальным решением квадратичной задачи программирования Equation 1. Однако, если любой компонент _λk отрицателен, а компонент не соответствует ограничению равенства, то минимизация не завершена. Алгоритм удаляет элемент, соответствующий самому отрицательному умножителю, и начинает новую итерацию.

Иногда, когда решатель заканчивает со всеми неотрицательными множителями Лагранжа, мера оптимальности первого порядка выше допуска, или ограничительный допуск не достигается. В этих случаях решатель пытается достичь лучшего решения, следуя процедуре перезапуска, описанной в [1]. В этой процедуре решатель отбрасывает текущий набор активных ограничений, ослабляет ограничения, накладываемые на бит, и создает новый набор активных ограничений, пытаясь решить проблему таким образом, чтобы избежать цикличности (неоднократно возвращаясь в то же состояние). При необходимости решатель может выполнить процедуру перезапуска несколько раз.

Примечание

Не пытайтесь остановить алгоритм раньше, задавая большие значения ConstraintTolerance и OptimalityTolerance опции. Обычно решатель выполняет итерацию, не проверяя эти значения, пока не достигнет потенциальной точки остановки, и только затем проверяет, выполняются ли допуски.

Иногда, active-set алгоритм может испытывать трудности с обнаружением, когда задача неограниченна. Эта проблема может возникнуть, если направление v несвязанности является направлением, где квадратичный член v'Hv = 0. Для числовой устойчивости и включения факторизации Холесского, active-set алгоритм добавляет небольшой, строго выпуклый член к квадратичной цели. Этот маленький термин заставляет целевую функцию быть ограниченной от - Inf. В этом случае active-set алгоритм достигает предела итерации вместо сообщения о том, что решение неограниченно. Другими словами, алгоритм останавливается с выходным флагом 0 вместо - 3.

Ссылки

[1] Гилл, П. Е., У. Мюррей, М. А. Сондерс и М. Х. Райт. Практическая процедура анти-циклирования для линейно ограниченной оптимизации. Math. Программирование 45 (1), август 1989, стр. 437-474.

Теплый старт

Когда вы запускаете quadprog или lsqlin 'active-set' алгоритм с теплым начальным объектом в качестве начальной точки, решатель пытается пропустить многие шаги Фазы 1 и Фазы 2. Объект теплого запуска содержит активный набор ограничений, и этот набор может быть правильным или близким к исправлению новой задачи. Поэтому решатель может избежать итераций, чтобы добавить ограничения к активному набору. Кроме того, начальная точка может быть близка к решению новой задачи. Для получения дополнительной информации см. optimwarmstart.

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