Изменения

Convex hull trick {{---}} один из методов оптимизации [[Динамическое_программирование | динамического программирования]], использующий идею [[Статические_выпуклые_оболочки:_Джарвис,_Грэхем,_Эндрю,_Чен,_QuickHull|выпуклой оболочки]]. Позволяет улучшить асимптотику решения некоторых задач, решемых методом динамического программирования, с <math>O(n^2)</math> до <tex>O(n\cdot\log(n))</tex>. Техника впервые появилась в 1995 году (задачу на нее предложили в USACO {{---}} национальной олимпиаде США по программированию). Массовую известность получила после IOI (международной олимпиады по программированию для школьников) 2002. ==Постановка примера Пример задачи, решаемой методом convex hull trick==Рассмотрим задачу на ДП:{{Задача|definition = Есть <math>n </math> деревьев с высотами a1<tex>a_1, a_2, a2\dots, … ana_n</tex> (в метрах). Требуется спилить их все, потратив минимальное количество монет на заправку бензопилы. Но пила устроена так, что она может спиливать только по <math>1</math> метру от дерева, к которому ее применили. Также после уменьшения высоты спиливаемого срубленного метра (любого дерева на 1 ее надо заправить) пилу нужно заправлять, платя за бензин определенное кол-во монет. Причем стоимость заправки бензина зависит от срубленных (полностью) деревьев. Если сейчас максимальный индекс срубленного дерева равен <tex>i</tex>, то цена заправки равна ci<tex>c_i</tex>. Изначально пила заправлена.И Также известны следующие ограничения : c[n] <tex>c_n = 0, a[1] a_1 = 1, a[i] a_i</tex> возрастают, c[i] <tex>c_i</tex> убывают.Изначально пила заправлена.(убывание и возрастание нестрогие)}}(Задача H с codeforcesСанкт-Петербургских сборов к РОИ 2016 <ref>[http://neerc.ifmo.ru/school/camp-2016/problems/20160318a.compdf Сайт с задачами Санкт-Петербургских сборов к РОИ 2016]</ref>)</noinclude><includeonly>{{#if: {{{neat|}}}|<div style="background-color: #fcfcfc; float:left;"><div style="background-color: #ddd;">'''Задача:'''</div><div style="border:1px dashed #2f6fab; padding: 8px; font-style: italic;">{{{definition}}}</div></div>|<table border="0" width="100%"><tr><td style="background-color: #ddd">'''Задача:'''</td></tr><tr><td style="border:1px dashed #2f6fab; padding: 8px; background-color: #fcfcfc; font-style: italic;">{{{definition}}}</td></tr></table>}}</includeonly> 

Сначала заметим важный факт : т.к. <tex>c[i]</tex> убывают (нестрого) и <tex>c[n] = 0</tex>, то все <tex>c[i]</tex> неотрицательны.Понятно, что нужно затратив минимальную стоимость срубить последнее (<tex>n</tex>-е) дерево, т.к. после него все деревья можно будет пилить рубить бесплатно (т.к. <tex>c[n] = 0</tex>). Посчитаем следующую динамику : <tex>dp[i] </tex> {{- --}} минимальная стоимость, заплатив которую будет срублено можно добиться того, что дерево номер <tex>i</tex> будет срублено. Тогда База динамики : <tex>dp[i1] = min(dp[j] + a[i] * c[j]) 0</tex>, т.к. изначально пила заправлена и высота первого дерева равна <math>1</math>, по всем j < iусловию задачи. То есть Переход динамики : понятно, что выгодно рубить сначала более дорогие и низкие деревья, а потом более высокие и дешвые дешевые (док-во этого факта оставляется читателям как несложное упражнение, т.к. эта идея относится скорее к теме [[Теорема_Радо-Эдмондса_(жадный_алгоритм)|жадных алгоритмов]], чем к теме данной статьи). Поэтому перед <tex>i</tex>-м деревом мы обязательно срубили какое-то <tex>j</tex>-е, причем <tex>j \leqslant i - 1</tex>. Тогда переберем Поэтому чтобы найти <tex>dp[i]</tex> нужно перебрать все <tex>1 \leqslant j \leqslant i - 1</tex> и попытаться использовать ответ для дерева номер <tex>j</tex>. Итак, пусть перед <tex>i</tex>-м деревом мы полностью срубили <tex>j < /tex>-е, причем высота <tex>i </tex>- го дерева составляет <tex>a[i]</tex>, а т.к. последнее дерево, которое мы срубили, имеет индекс предыдущего срубленного <tex>j</tex>, то стоимость каждого метра <tex>i</tex>-го дерева составит <tex>c[j]</tex>. Поэтому на сруб <tex>i</tex>-го деревамы потратим <tex>a[i] \cdot c[j]</tex> монет. Также не стоит забывать, что ситуацию, когда <tex>j</tex>-е дерево полностью срублено, мы получили не бесплатно, а за <tex>dp[j]</tex> монет.Итоговая формула пересчета : <tex>dp[i] = \min\limits_{j=1...i-1} (dp[j] + a[i] \cdot c[j])</tex>. Посмотрим на код вышеописанного решения: '''int''' <tex>\mathtt{simpleDP}</tex>('''int''' a[n], '''int''' c[n]) dp[1] = 0 dp[2] = dp[3] = .. . = dp[n] = <tex>\infty</tex> '''for''' i = 2 = 1..n dp[i] = <tex>+\infty</tex> '''for''' j = 1..i - 1 '''if''' (dp[j] + a[i] <tex>\cdot</tex> c[j] < dp[i]) dp[i] = dp[j] + a[i] <tex>\cdot</tex> c[j] '''return''' dp[n]Нетрудно видеть, что такая динамика работает за <tex>O(n^2)</tex>. ==Ключевая идея оптимизации==Для начала сделаем замену обозначений. Давайте обозначим <tex>dp[j]</tex> за <tex>b[j]</tex>, <tex>a[i]</tex> за <tex>x[i]</tex>, а <tex>c[j]</tex> за <tex>k[j]</tex>. Теперь формула приняла вид <tex>dp[i] = \min\limits_{j=0...i-1}(k[j] \cdot x[i] + b[j])</tex>. Выражение <tex>k[j] \cdot x + b[j]</tex> {{---}} это в точности уравнение прямой вида <tex>y = kx + b</tex>. Сопоставим каждому <tex>j</tex>, обработанному ранее, прямую <tex>y[j](x) = k[j] \cdot x + b[j]</tex>. Из условия «<tex>c[i]</tex> убывают <tex>\Leftrightarrow k[j]</tex> уменьшаются с номером <tex>j</tex>» следует то, что прямые, полученные ранее отсортированы в порядке убывания углового коэффициент. Давайте нарисуем несколько таких прямых : [[Файл:picture1convexhull.png]] Выделим множество точек <tex>(x_0, y_0)</tex> , таких что все они принадлежат одной из прямых и при этом нету ни одной прямой <tex>y’(x)</tex>, такой что <tex>y’(x_0) < y_0</tex>. Иными словами возьмем «выпуклую (вверх) оболочку» нашего множества прямых (её еще называют нижней огибающей множества прямых на плоскости). Назовем ее «<tex>y = convex(x)</tex>». Видно, что множество точек <math>(x, convex(x))</math> представляет собой выпуклую вверх функцию. ==Цель нижней огибающей множества прямых==Пусть мы его срубили отптимальным считаем динамику для <tex>i</tex>-го дерева. Его задает <tex>x[i]</tex>. Итак, нам нужно для данного <tex>x[i]</tex> найти <tex>\min\limits_{j=0..i-1}(k[j] \cdot x[i] + b[j]) = \min\limits_{j=0..i-1}(y[j](x[i]))</tex>. Это выражение есть <math>convex(x[i])</math>. Из монотонности угловых коэффицентов отрезков, задающих выпуклую оболочку, и их расположения по координатам x следует то, что отрезок, который пересекает прямую <tex>x = x[i]</tex>, можно найти [[Целочисленный_двоичный_поиск|бинарным поиском]]. Это потребует <tex>O(\log(n))</tex> времени на поиск такого <tex>j</tex>, что <tex>dp[i] = k[j] \cdot x[i] + b[j]</tex>. Теперь осталось научиться поддерживать множество прямых и быстро добавлять <tex>i</tex>-ю прямую после того, как мы посчитали <tex>b[i] = dp[i]</tex>. Воспользуемся идеей алгоритма построения выпуклой оболочки множества точек. Заведем 2 стека <tex>k[]</tex> и <tex>b[]</tex>, которые задают прямые в отсортированном порядке их угловыми коэффицентами и свободными членами. Рассмотрим ситуацию, когда мы хотим добавить новую (<tex>i</tex>-тую) прямую в множество. Пусть сейчас в множестве лежит <tex>sz</tex> прямых (нумерация с 1). Пусть <tex>(x_L, y_L)</tex> {{---}} точка пересечения <tex>sz - 1</tex>-й прямой множества и <tex>sz</tex>-й, а <tex>(x_R, y_R)</tex> {{---}} точка пересечения новой прямой, которую мы хотим добавить в конец множества и <tex>sz</tex>-й. Нас будут интересовать только их <tex>x</tex>-овые координаты <tex>x_L</tex> и <tex>x_R</tex>, соответственно. Если оказалось, что новая прямая пересекает <tex>sz</tex>-ю прямую выпуклой оболочки позже, чем <tex>sz</tex>-я <tex>sz - 1</tex>-ю, т.е. <tex>(x_L \geqslant x_R)</tex>, то <tex>sz</tex>-ю удалим из нашего множества, иначе - остановимся. Так будем делать, пока либо число прямых в стеке не станет равным 2, либо <tex>x_L</tex> не станет меньше <tex>x_R.</tex> Асимптотика : аналогично обычному алгоритму построения выпуклой оболочки, каждая прямая ровно <math>1</math> раз добавится в стек и максимум <math>1</math> раз удалится. Значит время работы перестройки выпуклой оболочки займет <tex>O(n)</tex> суммарно. [[Файл:picture2convexhull.png]][[Файл:picture3convexhull.png]] {{Теорема|id=th1239. |statement=Алгоритм построения нижней огибающей множества прямых корректен.|proof=Достаточно показать, что последнюю прямую нужно удалить из множества <tex>\Leftrightarrow</tex>, когда наша новая прямая пересекает ее в смысле денегточке с координатой по оси X, меньшей, чем последняя - предпоследнюю. Пусть <tex>Y(x) = Kx + B</tex> {{---}} уравнение новой прямой, <tex>y[i](x) способом= K[i]x + B[i]</tex> {{---}} уравнения прямых множества. Тогда просто т.к. <tex>K < K[sz]</tex>, то при <tex>x \in [- \infty; x_R] : y[sz](x) <= Y(x)</tex>, а т.к. <tex> K[sz] < K[sz - 1]</tex>, то при <tex>x \in [x_L; + \infty] : y[sz - 1](x) \geqslant y[sz](x)</tex>. Если <tex>x_L < x_R</tex>, то при <tex>x \in [x_L; x_R] : y[sz - 1] \geqslant y[sz](x) и Y(x) \geqslant y[sz](x)</tex>, т.е. на отрезке <tex>[x_L; x_R]</tex> прямая номер sz лежит ниже остальных и её нужно оставить в множестве. Если же <tex>x_L > x_R</tex>, то она ниже всех на отрезке <tex>[x_L; x_R] = \varnothing </tex>, т.е. её можно удалить из множества.}} ==Детали реализации:==Будем хранить 2 массива : <tex>front[]</tex> {{---}} <tex>x</tex>-координаты, начиная с которых прямые совпадают с выпуклой оболочкой (т.е. i-я прямая совпадает с выпуклой оболочкой текущего множества прямых при <tex>x</tex> <tex>\in</tex> <tex>[front[i]; front[i + 1])</tex> ) и <tex>st[]</tex> {{---}} номера деревьев, соответствующих прямым (т.е. <tex>i</tex>-я прямая множества, где <tex>i</tex> <tex>\in</tex> <tex>[1; sz]</tex> соответствует дереву номер <tex>sz[i]</tex>). Также воспользуемся тем, что <tex>x[i] = a[i] раз уменьшим высоту дерева </tex> возрастают (по условию задачи), а значит мы можем искать первое такое <tex>j</tex>, что <tex>x[i ] \geqslant front[j]</tex> не бинарным поиском, а методом двух указателей за <tex>O(n)</tex> операций суммарно. Также массив <math>front[]</math> можно хранить в целых числах, округляя х-координаты в сторону лежащих правее по оси x до ближайшего целого (*), т.к. на самом деле мы, считая динамику, подставляем в уравнения прямых только целые <tex>x[i]</tex>, а значит если <tex>k</tex>-я прямая пересекается с <tex>k+1</tex>-й в точке <tex>z +</tex> <tex>\alpha</tex> (<math>z</math>-целое, <tex>\alpha</tex> <tex>\in</tex> <tex>[0;1)</tex>), то мы будем подставлять в их уравнения <tex>z</tex> или <tex>z + 1</tex>. Каждый такой раз будем платить Поэтому можно считать, что новая прямая начинает совпадать с выпуклой оболочкой, начиная с <tex>x = z+1</tex> ==Реализация== '''int''' <tex>\mathtt{ConvexHullTrick}</tex>('''int''' a[n], '''int''' c[n]) st[1] = 1 front[1] = -<tex>\infty</tex><font color=green>// первая прямая покрывает все x-ы, начиная с -∞ </font> sz = 1 <font color=green>// текущий размер выпуклой оболочки </font> pos = 1 <font color=green>// текущая позиция первого такого j, что x[i] \geqslant front[st[j]] </font > '''for''' i = 2..n '''while''' (front[pos] < x[i]) <font color=green>// метод 1 указателя (ищем первое pos, такое что x[i] покрывается "областью действия" st[pos]-той прямой </font > pos = pos + 1 j = st[pos] dp[i] = K[j] <math>\cdot</math> a[i] + B[j] '''if''' (i < n) <font color=green>// если у нас добавляется НЕ последняя прямая, то придется пересчитать выпуклую оболочку </font > K[i] = c[i] <font color=green>// наши переобозначения переменных </font > B[i] = dp[i] <font color=green>// наши переобозначения переменных </font > x = -<tex>\infty</tex> '''while''' ''true'' j = st[sz] x = divide(B[j] - B[i], K[i] - K[j] ) <font color=green>// x-координата пересечения с последней прямой оболочки, округленное в нужную сторону (*) </font > '''if''' (x > from[sz]) '''break''' <font color=green>// перестаем удалять последнюю прямую из множества, если новая прямая пересекает ее позже, чем начинается ее "область действия" </font > sz = sz - 1<font color=green>// удаляем последнюю прямую, если она лишняя </font > st[sz + 1] = i front[sz + 1] = x <font color=green>// добавили новую прямую </font > sz = sz + 1 '''return''' dp[n] Здесь функция <tex>\mathtt{divide(a, b)}</tex> возвращает нужное(*) округление <tex>\frac{a}{b}</tex>. Приведем её код : '''int''' <tex>\mathtt{divide}</tex>('''int''' a, '''int''' b) delta = 0 '''if''' (a '''mod''' b ≠ 0) delta = 1 '''if''' ((a > 0 '''and''' b > 0) '''or''' (a < 0 '''and''' b < 0)) '''return''' [a / b] + delta '''return''' -[|a| / |b|]  Такая реализация будет работать за <math>O(n)</math>. ==Динамический convex hull trick== Заметим, что условия на возрастание/убывание <tex>k[i]</tex> на убывание/возрастание и <tex>x[i]</tex> выглядят достаточно редкими для большинства задач. Пусть в задаче таких ограничений нет. Первый способ борьбы с этой проблемой - отсортировать входные данные нужным образом, не испортив свойств задачи (пример : задача G c Санкт-Петербургских сборов к РОИ 2016 <ref>[http://neerc.ifmo.ru/school/camp-2016/problems/20160318a.pdf Сайт с задачами Санкт-Петербургских сборов к РОИ 2016]</ref>).  Но рассмотрим общий случай. По-прежнему у нас есть выпуклая оболочка прямых, с помощью которой мы за <tex>O(\log(n))</tex> можем найти <tex>dp[i]</tex>, но теперь вставку <tex>i</tex>-й прямой в оболочку уже нельзя выполнить описанным ранее способом за последующую заправку пилы<tex>O(1)</tex> в среднем. У нас есть выпуклая оболочка, наша прямая пересекает ее, возможно, «отсекая» несколько отрезков выпуклой оболочки в середине (рис. 4 : красная прямая - та, которую мы хотим вставить в наше множество). Более формально : теперь наша новая прямая будет ниже остальных при <tex>x \in [x_1; x_2]</tex>, где <tex>x_1, x_2 \in R</tex> - точки пересечения с некоторыми прямыми, причем <tex>x_2</tex> не обязательно равно <tex>+ \infty</tex>[[Файл:picture4convexhull.png]] Чтобы уметь вставлять прямую в множество будем хранить [[Красно-черное_дерево|двоичное дерево поиска]], в вершинах которого будут пары <tex>(k, st)</tex> = (коэффицент прямой, ее номер в глобальной нумерации). Когда приходит новая прямая, ищем последнюю прямую с меньшим угловым коэффицентом, чем у той прямой, которую мы хотим добавить в множество. Поиск такой прямой занимает <tex>O(\log(n))</tex>. Начиная с найденной прямой выполняем "старый" алгоритм (удаляем, пока текущая прямая множества бесполезна). ИтакИ симметричный алгоритм применяем ко всем прямым справа от нашей (удаляем правого соседа нашей прямой, пока она пересекает нас позже, чем своего правого соседа). Асимптотика решения составит <tex>O(\log(n))</tex> на каждый из <tex>n</tex> запросов «добавить прямую» + <tex>O(n\cdot\log(n))</tex> суммарно на сруб удаление прямых, т.к. по-прежнему каждая прямая не более одного раза удалится из множества, а каждое удаление из std::set занимает <tex>O(\log(n))</tex> времени. Итого <math>O(n\cdot\log(n))</math>. == Альтернативный подход ==Другой способ интерпретировать выражение <tex>dp[i] = \min\limits_{j=0...i-го дерева 1}(c[j] \cdot a[i] + dp[j])</tex> заключается в том, что мы заплатили будем пытаться свести задачу к стандартной выпуклой оболочке множества точек. Перепишем выражение средующим образом : <tex>dp[j] + a[i] \cdot c[j] = (dp[j], c[j]) \cdot (1, a[i]*)</tex>, т.е. запишем как скалярное произведение векторов <tex>v[j] = (dp[j], c[j])</tex> и <tex >u[i] = (1, a[i])</tex >. Вектора <tex >v[j] = (dp[j], c[j])</tex> хотелось бы организовать так, чтобы за <tex >O(\log(n))</tex> находить вектор, максимизирующий выражение <tex>v[j] \cdot u[i]</tex>. Посмотрим на рис. 
Нетрудно видеть5. Заметим интуитивно очевидный факт : красная точка (вектор) <tex>j</tex> не может давать более оптимальное значение <tex>v[j] \cdot u[i]</tex> одновременно чем обе синие точки. По этой причине нам достаточно оставить выпуклую оболочку векторов <tex>v[j]</tex>, а ответ на запрос {{---}} это поиск <tex>v[j]</tex>, максимизирующего проекцию на <tex>u[i]</tex>. Это задача поиска ближайшей точки выпуклого многоугольника (составленного из точек выпуклой оболочки) к заданной прямой (из <tex>(0, 0)</tex> в <tex>(1, что такая динамика работает a[i])</tex>). Ее можно решить за <tex>O(\log(n))</tex> двумя бинарными или одним тернарным поискомАсимптотика алгоритма по-прежнему составит <tex>O(n^2\cdot \log(n))</tex> [[Файл:picture5convexhull.png]] Докажем то, что описанный выше алгоритм корректен. Для этого достаточно показать, что если имеются <math>3</math> вектора <math>a, b, c</math>, расположенные как на рис. 5, т.е. точка <math>b</math> не лежит на выпуклой оболочке векторов <tex>0, a, b, c </tex> : <tex> \Leftrightarrow [a-b, b-c] < 0 </tex>, то либо <tex>(a, u[i])</tex> оптимальнее, чем <tex>(b, u[i])</tex>, либо <tex>(c, u[i])</tex> оптимальнее, чем <tex>(b, u[i])</tex>.{{Теорема|id=th12392. |statement=Если есть <tex>3</tex> вектора <tex>a, b, c</tex>, таких что <tex>[a-b, b-c] < 0</tex> то либо <math>(a, u) < (b, u)</math>, либо <math>(c, u) < (b, u)</math>, где вектор <math>u = (1; k)</math>.|proof=По условию теоремы известно, что <tex>[a-b, b-c] < 0 \Leftrightarrow (a_{x} - b_{x})\cdot(b_{y} - c_{y}) < (a_{y} - b_{y}) \cdot (b_{x} - c_{x})</tex> (*). Предположим (от противного), что <tex>(b, u) < (a, u) \Leftrightarrow b_{x} + k \cdot b_{y} < a_{x} + k \cdot a_{y} \Leftrightarrow (b_{x} - a_{x}) < k \cdot (a_{y} - b_{y})</tex> и при этом <tex>(b, u) < (c, u) \Leftrightarrow b_{x} + k \cdot b_{y} < c_{x} + k \cdot c_{y} \Leftrightarrow (c_{x} - b_{x}) > k \cdot (b_{y} - c_{y})</tex>. Подставим эти неравенства в (*). Получим цепочку неравенств : <tex>k \cdot (a_{y} - b_{y})</tex><tex> \cdot (c_{y} - b_{y}) = k</tex><tex> \cdot (b_{y} - a_{y}) \cdot </tex><tex>(b_{y} - c_{y})</tex> <tex> < (a_{x} - b_{x})</tex><tex> \cdot (b_{y} - c_{y})</tex><tex> < (a_{y} - b_{y}) \cdot </tex><tex>(b_{x} - c_{x})</tex> <tex>< k \cdot (a_{y} - b_{y})</tex><tex> \cdot (c_{y} - b_{y})</tex>. Получили противоречие : <tex>k \cdot (a_{y} - b_{y}) \cdot (c_{y} - b_{y}) < k \cdot (a_{y} - b_{y}) \cdot (c_{y} - b_{y})</tex>. Значит предположение неверно, чтд.}} Из доказанной теоремы и следует корректность алгоритма. ==См. также== *[[:Статические_выпуклые_оболочки:_Джарвис,_Грэхем,_Эндрю,_Чен,_QuickHull|Выпуклая оболочка]] *[[:Динамическое_программирование|Динамическое программирование]] == Примечания ==<references/> [[Категория:Дискретная математика и алгоритмы]][[Категория: Динамическое программирование]][[Категория: Способы оптимизации методов динамического программирования]]