Convex hull trick
Содержание
Что такое convex hull trick
Convex hull - выпуклая оболочка; Convex hull trick - один из методов оптимизации динамического программирования, использующий идею выпуклой оболочки. Позволяет улучшить ассимптотику решения некоторых задач, решемых методом динамического программирования с
до . Техника впервые появилась в 1995 году (задачу на нее предложили в USACO - национальной олимпиаде США по программированию). Массовую известность получила после IOI (международной олимпиады по программированию для школьников) 2002Пример задачи, решаемой методом convex hull trick
Рассмотрим задачу на ДП:
Есть n деревьев с высотами(в метрах). Требуется спилить их все, потратив минимальное количество монет на заправку бензопилы. Но пила устроена так, что она может спиливать только по 1 метру от дерева, к которому ее применили. Также после срубленного метра (любого дерева) пилу нужно заправлять, платя за бензин определенной кол-во монет. Причем стоимость бензина зависит от срубленных (полностью) деревьев. Если сейчас максимальный индекс срубленного дерева равен i, то цена заправки равна ci. Изначально пила заправлена. Также известны следующие ограничения : возрастают, убывают. Изначально пила заправлена. (убывание и возрастание нестрогие)
(Задача H отсюда : http://neerc.ifmo.ru/school/camp-2016/problems/20160318a.pdf)
Наивное решение
Сначала заметим важный факт : т.к. c[i] убывают(нестрого) и c[n] = 0, то все c[i] неотрицательны. Понятно, что нужно затратив минимальную стоимость срубить последнее (
-е) дерево, т.к. после него все деревья можно будет рубить бесплатно (т.к. ). Посчитаем следующую динамику : - минимальная стоимость, заплатив которую можно добиться того, что дерево номер будет срублено. База динамики : , т.к. изначально пила заправлена и высота первого дерева равна 1, по условию задачи. Переход динамики : понятно, что выгодно рубить сначала более дорогие и низкие деревья, а потом более высокие и дешевые (док-во этого факта оставляется читателям как несложное упражнение, т.к. эта идея относится скорее к теме жадных алгоритмнов, чем к теме данной статьи). Поэтому перед -м деревом мы обязательно срубили какое-то -е, причем . Поэтому чтобы найти нужно перебрать все и попытаться использовать ответ для дерева намер . Итак, пусть перед -м деревом мы полностью срубили -е, причем высота -го дерева составляет , а т.к. последнее дерево, которое мы срубили имеет индекс , то стоимость каждого метра -го дерева составит . Поэтому на сруб -го дерева мы потратим монет. Также не стоит забывать, ситуацию, когда -е дерево полностью срублено, мы получили не бесплатно, а за монет. Итогвая формула пересчета : . Посмотрим на код выше описанного решения:dp[1] = 0 dp[2] = dp[3] = ... = dp[n] =for i = 1..n-1 { dp[i] = + for j = 0..i-1 { if (dp[j] + a[i] * c[j] < dp[i]) dp[i] = dp[j] + a[i] * c[j] } }
Нетрудно видеть, что такая динамика работает за
.Ключевая идея оптимизации
1) Сделаем замену обозначений. Давайте обозначим
за , за , а за .2) Теперь формула приняла вид
. Выражение - это в точности уравнение прямой вида .3) Сопоставим каждому
, обработанному ранее, прямую . Из условия « убывают уменьшаются с номером » следует то, что прямые, полученные ранее отсортированы в порядке убывания углового коэффициент. Давайте нарисуем несколько таких прямых :4) Выделим множество точек
, таких что все они принадлежат одной из прямых и при этом нету ни одной прямой , такой что . Иными словами возьмем «выпуклую (вверх) оболочку» нашего множества прямых (её еще называют нижней ошибающей множества прямых на плоскости). Назовем ее « ». Видно, что множество точек (x, convex(x)) представляет собой выпуклую вверх функцию.Для чего нужна нижняя ошибающая множества прямых
Пусть мы считаем динамику для
-го дерева. Его задает . Итак, нам нужно для данного найти . Это выражение есть convex(x[i]). Из монотонности угловых коэффицентов отрезков, задающих выпуклую оболочку, и их расположения по координаты x следует то, что отрезок, который пересекает прямую , можно найти бинарным поиском. Это потребует времени на поиск такого j, что dp[i] = k[j] * x[i] + b[j]. Теперь осталось научиться поддерживать множество прямых и быстро добавлять -ю прямую после того, как мы посчитали .Воспользуемся идеей алгоритма построения выпуклой оболочки множества точек. Заведем 2 стека
и , которые задают прямые в отсортированном порядке их угловыми коэффицентами и свободными членами. Рассмотрим ситуацию, когда мы хотим добавить новую ( -тую) прямую в множество. Найдем точки пересечения с последними 2мя прямыми из стека. Нас будут интересовать только их -овые координаты. Назовем их и , соответственно. Если оказалось, что новая прямая пересекает предпоследнюю прямую выпуклой оболочки позже, чем последнюю , то последнюю удалим из нашего множества. Так будем делать, пока либо кол-во прямых в стеке не станет равным 2, либо не станет меньшеАсимптотика : аналогично обычному алгоритму построения выпуклой оболочки, каждая прямая ровно 1 раз добавится в стек и максимум 1 раз удалится. Значит время работы перестройки выпуклой оболочки займет
суммарно.Корректность : достаточно показать, что последнюю прямую нужно удалить из множества т.и т.т., когда она наша новая прямая пересекает ее в точке с координатой по оси X, меньшей, чем предпоследнюю.
План доказательства : пусть новая прямая пересекает последнюю прямую множества позже, чем предпоследнюю (рис.2 - красная прямая новая, фиолетовая - предпоследняя, желтая - последняя), тогда найдется такой отрезок
, что последняя (желтая) прямая при лежит ниже всех остальных и ее следует оставить в множестве. А если новая прямая пересекает последнюю прямую множества раньше, чем предпоследнюю (рис.1), последняя прямая при любых x лежит выше какой-то другой прямой множества и значит ее нужно удалить.Более формально. Пусть
Детали реализации:
Будем хранить 2 массива (имитирующих стеки) :
и - начало (по x) соответствующей прямой выпуклой оболочки и номер этой прямой (в глобальной нумерации). Также воспользуемся тем, что возрастают (по условию задачи), а значит мы можем искать первое такое , что не бинарным поиском, а методом 2х указателей за суммарно. Также массив front[] можно хранить в целых числах, округляя х-координаты в сторону лежащих правее по оси x до ближайшего целого (*). Почему так? На самом деле мы, считая динамику, подставляем в уравнения прямых только целые , а значит если -я прямая пересекается с -й в точке (z-целое, ), то мы будем подставлять или . Поэтому можно считать, что новая прямая имеет "область действия", начиная сРеализация
void Convex-hull-trick st[0] = 0 from[0] = -// первая прямая покрывает все x-ы, начиная с -∞ sz = 1 // текущий размер выпуклой оболочки pos = 0 // текущая позиция первого такого j, что x[i] >= front[st[j]] for i = 1..n-1 { while (front[pos] < x[i]) // метод 1 указателя (ищем первое pos, такое что x[i] покрывается "областью действия" st[pos]-той прямой ++pos j = st[pos] dp[i] = K[j] * a[i] + B[j] if (i < n - 1) { // если у нас добавляется НЕ последняя прямая, то придется пересчитать выпуклую оболочку K[i] = c[i] // наши переобозначения переменных B[i] = dp[i] // наши переобозначения переменных x = - while (true) { j = st[sz - 1] x = divide(B[j] - B[i], K[i] - K[j]) // x-координата пересечения с последней прямой оболочки, округленное в нужную сторону (*) if (x > from[sz - 1]) break // перестаем удалять последнюю прямую из множества, если новая прямая пересекает ее позже, чем начинается ее "область действия" --sz // удаляем последнюю прямую, если она лишняя } st[sz] = i from[sz++] = x // добавили новую прямую } }
(Здесь функция divide(a, b) возвращает нужное(*) округление a / b) Такая реализация будет работать за O(n).
Динамический convex hull trick
Заметим, что условия на прямые, что http://neerc.ifmo.ru/school/camp-2016/problems/20160318a.pdf).
Но рассмотрим общий случай. Наша задача поменялась следующим образом : по-прежнему у нас есть выпуклая оболочка, имея которую мы за или быстрее можем найти , но теперь вставку i-й прямой в оболочку уже нельзя выполнить старым способом за (в среднем). У нас есть выпуклая оболочка, наша прямая пересекает ее, возможно, «отрезая» несколько отрезков выпуклой оболочки в середине (рис. 4).
Т.е. нужно уметь быстро (за
?) назодить, после какой прямой стоит пытаться вставить текущую (красную рис.4) прямую и удалять лишние справа, начиная с нее, потом проводить аналогичные операции слева. Итак, давайте хранить (или любой аналог в других языках) пар = <коэффицент прямой, ее номер в глобальной нумерации>. Когда приходит новая прямая, делаем lower_bound - 1 в сете, т.е. ищем ближайшую прямую с меньшим углом наклона, и начиная с нее повторяем старый алгоритм (удаляем, пока прямая бесполезная). И симметричный алгоритм применяем ко всем прямым справа от нашей. Асимптотика решения составит на каждый из n запросов «добавить прямую» + суммарно на удаление прямых. Итого .Альтернативный подход
Другой способ интерпретировать выражение
заключается в следующем: давайте перепишем выражение , т.е. запишем как скалярное произведение векторов и . Вектора хотелось бы организовать так, чтобы за находить максимизирующий выражение . Посмотрим на рис. 5. Заметим довольно очевидный факт : красная точка(вектор) не может давать более оптимальное значение одновременно чем обе синие точки, т.к. - это на самом деле проекция вектора на . По этой причине нам достаточно оставить выпуклую оболочку векторов , а ответ на запрос - это поиск , максимизирующего проекцию на . Это задача поиска ближайшей точки выпуклого многоугольника (составленного из точек выпуклой оболочки) к заданной прямой (из в ). Ее можно решить за O(logn) двумя бинарными или одним тернарным поиском Асимптотика алгоритма по-прежнему составит