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

    Я долго готовился и собирал материал, надеюсь в этот раз получилось лучше. Эту статью посвящаю основным методам решения задач математической оптимизации с ограничениями, так что если вы слышали, что симплекс-метод — это какой-то очень важный метод, но до сих пор не знаете, что он делает, то возможно эта статья вам поможет.

    P. S. Статья содержит математические формулы, добавленные макросами хабраредактора. Говорят, что они иногда не отображаются. Также есть много анимаций в формате gif.

    Преамбула


    Задача математической оптимизации — это задача вида «Найти в множестве $\mathcal{K}$ элемент $x^*$ такой, что для всех $x$ из $\mathcal{K}$ выполняется $f(x^*)\leq f(x)$», что в научной литературе скорее будет записано как-то так

    $\begin{array}{rl} \mbox{минимизировать } & f(x) \\ \mbox{при условии } & x\in \mathcal{K}. \end{array} $


    Исторически так сложилось, что популярные методы такие как градиентный спуск или метод Ньютона работают только в линейных пространствах (причем желательно простых, например $\mathbb{R}^n$). На практике же часто встречаются задачи, где нужно найти минимум не в линейном пространстве. Например нужно найти минимум некоторой функции на таких векторах $x=(x_1, \ldots, x_n)$, для которых $x_i\geq 0$, это может быть обусловлено тем, что $x_i$ обозначают длины каких-либо объектов. Или же например, если $x$ представляют координаты точки, которая должна быть на расстоянии не больше $r$ от $y$, т. е. $\|x-y\|\leq r$. Для таких задач градиентный спуск или метод Ньютона уже напрямую не применить. Оказалось, что очень большой класс задач оптимизации удобно покрывается «ограничениям», подобными тем, что я описал выше. Иначе говоря, удобно представлять множество $\mathcal{K}$ в виде системы равенств и неравенств

    $ \begin{array}{l} g_i(x)\leq 0,~1\leq i \leq m, \\ h_i(x)=0,~1\leq i\leq k. \end{array} $


    Задачи минимизации над пространством вида $\mathbb{R}^n$ таким образом стали условно называть «задачами без ограничений» (unconstrained problem), а задачи над множествами, заданными наборами равенств и неравенств — «задачами с ограничениями» (constrained problem).

    Технически, совершенно любое множество $\mathcal{K}$ можно представить в виде одного равенства или неравенство с помощью индикатор-функции, которая определяется как

    $I_\mathcal{K}(x)= \begin{cases} 0, & x\notin \mathcal{K} \\ 1, & x\in \mathcal{K}, \end{cases} $


    однако такая функция не обладает разными полезными свойствами (выпуклость, дифференцируемость и т. п.). Тем не менее, часто можно представить $\mathcal{K}$ в виде нескольких равенств и неравенств, каждое из которых такими свойствами обладает. Основная теория подведена под случай

    $\begin{array}{rl} \mbox{минимизировать } & f(x) \\ \mbox{при условии } & g_i(x)\leq 0,~1\leq i\leq m \\ & Ax=b, \end{array} $


    где $f, g_i$ — выпуклые (но не обязательно дифференцируемые) функции, $A$ — матрица. Для демонстрации работы методов я буду использовать два примера:

    1. Задача линейного программирования

      $$display$$\begin{array}{rl} \mbox{минимизировать } & -2&x~~~- &y \\ \mbox{при условии } &-1.0 & ~x -0.1 & ~y \leq -1.0 \\ & -1.0 & ~x + ~0.6 &~ y \leq -1.0 \\ & -0.2 & ~x + ~1.5 & ~y \leq -0.2\\ & ~0.7 &~x + ~0.7 & ~y \leq 0.7\\ &~2.0 & ~x -0.2 & ~y \leq 2.0\\ &~0.5 & ~x -1.0 & ~y \leq 0.5\\ &-1.0 & ~x -1.5& ~ y \leq -1.0\\ \end{array} $$display$$


      По сути эта задача состоит в том, чтобы найти самую дальнюю точку многоугольника в направлении (2, 1), решение задачи — точка (4.7, 3.5) — самая «правая» в многоугольнике). А вот собственно и сам многоугольник

    2. Минимизация квадратичной функции с одним квадратичным ограничением

      $\begin{array}{rl} \mbox{минимизировать } & 0.7 (x - y)^2 + 0.1 (x + y)^2\\ \mbox{при условии } & (x-4)^2+(y-6)^2\leq 9 \end{array} $




    Симплекс-метод


    Из всех методов, которые я покрываю этим обзором, симплекс-метод наверно является самым известным. Метод был разработан специально для линейного программирования и единственный из представленных достигает точного решения за конечное число шагов (при условии, что для вычислений используется точная арифметика, на практике это обычно не так, но в теории возможно). Идея симплекс-метода состоит из двух частей:

    1. Системы линейных неравенств и равенств задают многомерные выпуклые многогранники (политопы). В одномерном случае это точка, луч или отрезок, в двумерном — выпуклый многоугольник, в трехмерном — выпуклый многогранник. Минимизация линейной функции — это по сути нахождение самой «дальней» точки в определенном направлении. Думаю, интуиция должна подсказывать, что этой самой дальней точкой должна быть некая вершина, и это действительно так. В общем случае, для системы из $m$ неравенств в $n$-мерном пространстве вершина — это точка, удовлетворяющая системе, для которой ровно $n$ из этих неравенств обращаются в равенства (при условии, что среди неравенств нет эквивалентных). Таких точек всегда конечное число, хоть их и может быть очень много.
    2. Теперь у нас есть конечный набор точек, вообще говоря можно их просто взять и перебрать, то есть сделать что-то такое: для каждого подмножества из $n$ неравенств решить систему линейных уравнений, составленных на выбранных неравенствах, проверить, что решение подходит в исходную систему неравенств и сравнить с другими такими точками. Это довольно простой неэффективный, но рабочий метод. Симплекс-метод вместо перебора двигается от вершины к вершине по ребрам таким образом, чтобы значений целевой функции улучшалось. Оказывается, если у вершины нет «соседей», в которых значений функции лучше, то она оптимальна.

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

    $\begin{array}{rl} \mbox{минимизировать } & s \\ \mbox{при условии } & g_i(x)\leq s,~1\leq i\leq m \\ \end{array} $


    Если для решения этой задачи $x^*, s^*$ такое, что $s^*\leq 0$, то выполняется $g_i(x^*)\leq s\leq 0$, иначе исходная задача вообще задана на пустом множестве. Чтобы решить вспомогательную задачу, можно также использовать симплекс-метод, начальной же точкой можно взять $s=\max_ig_i(x)$ с произвольным $x$. Нахождение начальной точки можно условно назвать первой фазой метода, нахождение решение исходной задачи можно условно назвать второй фазой метода.

    Траектория двухфазового симплекс-метода
    Траектория была сгенерирована с помощью scipy.optimize.linprog.


    Проективный градиентный спуск


    Про градиентный спуск я недавно писал отдельную статью, в которой в том числе кратко описал и этот метод. Сейчас этот метод вполне себе живой, но изучается как часть более общего проксимального градиентного спуска. Сама идея метода совсем банальна: если мы применяем градиентный спуск к выпуклой функции $f$, то при правильном выборе параметров получаем глобальный минимум $f$. Если же после каждого шага градиентного спуска корректировать полученную точку, взяв вместо нее её проекцию на замкнутое выпуклое множество $\mathcal{K}$, то в результате мы получим минимум функции $f$ на $\mathcal{K}$. Ну или более формально, проективный градиентный спуск — это алгоритм, который последовательно вычисляет

    $ \begin{cases} x_{k+1} = y_k - \alpha_k\nabla f(y_k) \\ y_{k+1} = P_{\mathcal{K}}(x_{k+1}), \end{cases} $


    где

    $ P_{\mathcal{K}}(x)=\mbox{argmin}_{y\in\mathcal{K}}\|x-y\|. $


    Последнее равенство определяет стандартный оператор проекции на множество, по сути это функция, которая по точке $x$ вычисляет ближайшую к ней точку множества $\mathcal{K}$. Роль расстояния здесь играет $\|\ldots\|$, стоит отметить, что здесь можно использовать любую норму, тем не менее проекции с разными нормами могут отличаться!

    На практике проективный градиентный спуск используется только в особых случаях. Основная его проблема состоит в том, что вычисление проекции может быть еще более сложной задачей, чем исходная, а её нужно вычислять много раз. Самый распространенный случай, для которого удобно применять проективный градиентный спуск — это «коробочные ограничения», которые имеют вид

    $ \ell_i\leq x_i\leq r_i,~~1\leq i\leq n. $


    В этом случае проекция вычисляется очень просто, по каждой координате получается

    $ [P_{\mathcal{K}}(x)]_i= \begin{cases} r_i, & x_i>r_i \\ \ell_i, & x_i < \ell_i \\ x_i, & x_i\in [\ell_i, r_i]. \end{cases} $


    Применение проективного градиентного спуска для задач линейного программирования совершенно бессмысленно, тем не менее если это все-таки сделать, то выглядеть будет как-то так

    Траектория проективного градиентного спуска для задачи линейного программирования


    А вот как выглядит траектория проективного градиентного спуска для второй задачи, если

    выбирать большой размер шага


    и если

    выбирать небольшой размер шага


    Метод эллипсоидов


    Этот метод примечателен тем, что является первым полиномиальным алгоритмом для задач линейного программирования, его можно считать многомерным обобщением метода бисекции. Я начну с более общего метода разделяющей гиперплоскости:

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

    Для задач оптимизации построение «разделяющей гиперплоскости» основано на следующем неравенстве для выпуклых функций

    $ f(y)\geq f(x)+\nabla f(x)^T(y-x). $


    Если зафиксировать $x$, то для выпуклой функции $f$ полупространство $\nabla f(x)^T(y-x)\geq 0$ содержит только точки со значением не меньше, чем в точке $x$, а значит их можно отсечь, так как эти точки не лучше, чем та, что мы уже нашли. Для задач с ограничениями можно аналогичным образом избавиться от точек, которые гарантированно нарушают какое-то из ограничений.

    Самый простой вариант метода разделяющей гиперплоскости — это просто отсекать полупространства без добавления каких-либо точек. В результате на каждом шаге у нас будет некий многогранник. Проблема этого метода в том, что количество граней многогранника скорее всего будет возрастать от шага к шагу. Более того, оно может расти экспоненциально.

    Метод эллипсоидов собственно на каждом шаге хранит эллипсоид. Точнее, после проведения гиперплоскости строится эллипсоид минимального объема, который содержит одну из частей исходного. Этого получается добиться за счет добавления новых точек. Эллипсоид всегда можно задать положительно определенной матрицей и вектором (центром эллипсоида) следующим образом

    $ \mathcal{E}(P, x)=\{z~|~(z-x)^TP(z-x)\leq 1\}. $


    Построение минимального по объему эллипсоида, содержащего пересечение полупространства и другого эллипсоида, можно осуществить с помощью в меру громоздких формул. К сожалению на практике этот метод оказался все еще на так хорош, как симплекс-метод или метод внутренней точки.

    А вот собственно как он работает для

    линейного программирования


    и для

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


    Метод внутренней точки


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

    Базовая идея метода — замена ограничений на штраф в виде так называемой барьерной функции. Функция $F:Int\mathcal{K}\rightarrow \mathbb{R}$ называется барьерной функцией для множества $\mathcal{K}$, если

    $ F(x)\rightarrow +\infty~\mbox{при } x\rightarrow\partial \mathcal{K}. $


    Здесь $Int\mathcal{K}$ — внутренность $\mathcal{K}$, $\partial\mathcal{K}$ — граница $\mathcal{K}$. Вместо исходной задачи предлагается решать задачу

    $ \mbox{минимизировать по }x~~\varphi(x, t)=tf(x)+F(x). $


    $F$ и $\varphi$ заданы только на внутренности $\mathcal{K}$ (по сути отсюда и название), свойство барьера гарантирует, что у $\varphi$ минимум по $x$ существует. Более того, чем больше $t$, тем больше влияние $f$. При достаточно разумных условиях можно добиться того, что если устремить $t$ к бесконечности, то минимум $\varphi$ будет сходиться к решению исходной задачи.

    Если множество $\mathcal{K}$ задано в виде набора неравенств $g_i(x)\leq 0,~1\leq i\leq m$, то стандартным выбором барьерной функции является логарифмический барьер

    $ F(x)=-\sum_{i=1}^m\ln(-g_i(x)). $


    Точки минимума $x^*(t)$ функции $\varphi(x, t)$ для разных $t$ образует кривую, которую обычно называют центральный путь, метод внутренний точки как бы пытается следовать этому пути. Вот так он выглядит для

    Примера с линейным программированием

    Аналитический центр — это просто $x^*(0)$

    Наконец сам метод внутренней точки имеет следующий вид

    1. Выбрать начальное приближение $x_0$, $t_0>0$
    2. Выбрать новое приближение методом Ньютона

      $ x_{k+1}=x_k-[\nabla^2_x \varphi(x_k, t_k)]^{-1}\nabla_x \varphi(x_k, t_k) $


    3. Увеличить $t$

      $ t_{k+1}=\alpha t_k,~~\alpha>1 $



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

    Задача линейного программирования

    Прыгающая черная точка — это $x^*(t_k)$, т.е. точка, к которой мы пытаемся приблизиться шагом метода Ньютона на текущем шаге.

    Задача квадратичного программирования

    Share post

    Comments 20

      +1
      Жаль, что материал на такую тему исключительно справочный… Есть ли ссылки на адекватные программные реализации? (Без симплекса — его и так много)
        0
        Смотря для каких целей. Если для прикладных — не стоит задумываться о том, как и какой метод реализован, лучше просто использовать готовый пакет оптимизации типа гуроби, cvxopt или scipy.optimize. Есть например хороший интерфейс на питон — cvxpy. В matlab точно что-то есть, скорее всего и в подобных математических продуктах тоже.

        Если нужно «пощупать» руками, я могу дать свои реализации, по которым анимации строил.
          0
          Спасибо. Я много лет использую библиотеку PuLP, вдруг пригодится… Кстати, какой метод реализован — на больших размерностях бывает важно… Во всяком случае, если вы хотите знать сколько времени (хотя бы порядок) может занять сам расчет. Вообще, у меня есть небольшая претензия к названию статьи (излишне обобщили), но это мелочи с учётом того, что вы продемонстрировали ;) Спасибо
          0
          В MATLAB есть пакет CVX, который является удобной надстройкой над солверами. Он же портирован в R.
          0
          Спасибо за статью. Добавьте, пожалуйста, список используемой литературы.
            0
            Отличная статья!
            По методу внутренней точки для линейного программирования:
            — логарифмический барьер — правильный? (шаг метода Ньютона генерирует точку, которая остается внутри нашего множества)
            — шаг метода L-BFGS с логарифмическим барьером генерирует точку, которая остается внутри нашего множества?
            — какие обычно значения берутся для t0 и a?
              0
              Логарифмический барьер — правильный. К сожалению тут есть небольшие тонкости с выбором начальных данных. Готового рецепта на все случаи жизни я к сожалению не готов дать. Когда готовил анимации у меня возникла проблема, что Ньютон разваливался. Тогда я это решил методом научного тыка.

              Про BFGS и все его вариации не могу ничего сказать, лично я пробовал применять градиентный спуск вместо Ньютона — ничего хорошего не выходило. Всегда можно самому попробовать, но думаю, что уже пробовали не один раз. Метод внутренней точки придумали давно и Ньютон там засел намертво.
              0
              Что значит, Ньютон разваливался?
              Еще не пробовал. Ищу методы с адекватной сложностью. Ваша статья очень помогла мне.
              Какие значения t0 и a вы брали для демонстрации?
                0
                Я использовал t=0.1, alpha=2. В некоторых экспериментах Ньютон на первом шаге уходил за пределы области. Студенты, которым давал задание реализовать метод внутренней точки, говорил, что тоже иногда Ньютон не выдавал точку внутри. В статье правку сделал
                0
                А симплекс-метод работает только для линейных ограничений и линейных функций?
                  0
                  Да
                    0
                    Вставлю свои пять копеек: для непрерывных величин (тип float и им подобные) и некоторых случаев булевого линейного программирования
                    0
                    Если не ошибаюсь, проективные градиенты и подобные ему — первого порядка сходимости, поэтому на практике, c общими нелинейными функциями, работают медленнее, чем, например, SQP, где локально решается задача квадратичного программирования. Я сам это наблюдал, когда использовал его с коробочными ограничениями. Если не присобачить поправку второго порядка, сходиться будет как черепаха.
                    А вообще классная тема. Обогнали меня, сам хотел что-то такое в перспективе написать!:)
                      0
                      Если говорить о количестве итераций до достижения определенной точности, то да — чем выше порядок метода, тем быстрее он сходится. Однако есть и другая сторона медали — чем выше порядок метода, тем сложнее одна итерация, и тем сложнее применять эти методы в пространствах очень больших размерностей. Например, градиентный спуск очень распространен в нейронных сетях, где обычно не меньше миллиона параметров, причем в основном его «стохастический» вариант, который очень хорошо укладывается в парадигму ML. Применять его в пространствах размера < 10000 практически бессмысленно, есть достаточно методов, которые будут работать быстрее.
                      0
                      Если не секрет, не подскажете, каким софтом рисовали анимацию сходящихся итераций для этого поста? Очень-очень полезная вещь!
                        0
                        Не секрет, это python/matplotlib. Я хотел приложить ноутбук к статье, но в итоге недооформил адекватно. Сейчас могу только выложить черновик
                          0
                          Спасибо!
                        0
                        Уточните, пожалуйста, обоснованность «Выбрать новое приближение методом Ньютона». Вы его применяете так, будто хотите найти не минимум, а 0 от функции fi(x,t).
                          0
                          Обратите внимание, что метод Ньютона применяется для градиента fi(x, t), т.е. мы действительно ищем 0 градиента функции, что соответствует минимуму самой функции в случае, когда она выпукла и дифференцируема, у нас как раз такой.
                            +1
                            А, действительно, вижу. У вас dx = -grad(fi) / grad(grad(fi)), а не dx = -fi / grad(fi). Разобрался, спасибо

                        Only users with full accounts can post comments. Log in, please.