Троичный поиск

Под каким углом бросить теннисный мяч, чтобы он улетел дальше всего? Без воздуха ответ знает любой школьник — 45°. С воздухом формулы нет: дальность считают, моделируя полёт шаг за шагом, и одна такая проба — тысячи шагов. Зато известно главное: если поднимать угол от нуля, дальность сначала растёт, а потом падает, и вершина у кривой одна. Значит, искать можно так же, как двоичным поиском, отбрасывая куски отрезка. Выберите снаряд и способ и посмотрите, сколько бросков понадобится, чтобы найти угол с точностью до сотой градуса.

Теннисный мяч, брошенный со скоростью 30 м/с с высоты метра, дальше всего летит под углом 38,02° — на 40,93 метра. Троичный поиск находит этот угол за 46 полётов, золотое сечение — за 21. Волан из перьев, который тормозит гораздо сильнее, надо бросать ещё положе, под 24,39°. И даже ядру воздух немного мешает: под 42,01° оно летит на 13 сантиметров ближе, чем пролетело бы в пустоте под 42,09°. Ровно 45° не годятся ни для кого: ядро выпускают с высоты 2,1 метра, а мяч и волан тормозит воздух.

Идея

Функцию на отрезке называют унимодальной, если она до некоторой точки строго растёт, а после неё строго убывает. Дальность броска от угла именно такая. Ни формула, ни производная для поиска не нужны — только возможность вычислить функцию в любой точке.

Возьмём внутри отрезка две точки m₁ < m₂ и сравним значения. Если в m₂ больше, вершина не может лежать левее m₁: слева от вершины функция растёт, и тогда в m₁ было бы меньше, чем в m₂, только если обе точки на подъёме. Отбрасываем кусок левее m₁. Если больше в m₁ — отбрасываем кусок правее m₂. Отрезок сжимается, а вершина остаётся внутри.

Классический троичный поиск ставит точки на третях отрезка и за два вычисления отбрасывает треть. Отсюда и название: отрезок делят на три части. В сцене это видно по нижней полосе: две точки, красная черта под выброшенной частью и подсвеченный остаток.

Почему вершина не теряется

Обратите внимание на слово «строго». Оно понадобится, когда речь пойдёт о площадках: если функция где-то постоянна, равенство значений ничего не говорит, и теорема не работает.

Из теоремы следует и конечность. Каждый шаг оставляет 2/3 отрезка, после k шагов длина равна (2/3)^k·L, и до точности ε нужно ⌈log_1,5(L/ε)⌉ шагов. Для угла от 0 до 90° с точностью 0,01° это 23 шага, то есть 46 полётов.

Золотое сечение: одно вычисление на шаг

Троичный поиск расточителен. После шага одна из двух точек оказывается внутри нового отрезка, и её значение уже известно, но на следующем шаге её место — не на трети, и значение выбрасывают. А что если поставить точки так, чтобы старая точка оказалась ровно там, где нужна новая?

Пусть отрезок — [0, 1], а точки стоят симметрично: a = 1 − r и b = r. Отбросим правый кусок: новый отрезок — [0, r], и в нём точки должны стоять там же в том же масштабе, то есть в r·(1 − r) и r·r. Старая точка a = 1 − r подойдёт на роль новой правой, если

Это золотое сечение, отсюда и имя метода: поиск золотым сечением. Каждый шаг оставляет 0,618 отрезка, но стоит одного вычисления, а не двух. Включите в сцене «Золотое сечение»: зелёный кружок отмечает точку, значение которой взято с прошлого шага.

Честнее всего сравнивать способы по тому, во сколько раз сжимается отрезок на одно вычисление функции:

Чтобы сжать отрезок в десять раз, троичному поиску нужно 11,4 вычисления, дихотомии — 6,6, золотому сечению — 4,8. Разница в 2,4 раза не зависит ни от функции, ни от точности: она во всём, что занимает время. Есть и ещё более экономный вариант — поиск Фибоначчи: если число вычислений известно заранее, точки ставят в отношениях соседних чисел Фибоначчи, и последний шаг выходит немного удачнее. Золотое сечение — его предел, когда вычислений много.

Код

Листинг ищет угол для ядра, мяча и волана тремя способами, считает полёты и проверяет себя формулой для броска в пустоте: tg θ = v/√(v² + 2gh). Полёт — метод Рунге — Кутты с шагом 2 мс, сопротивление воздуха квадратичное. Каждую функцию вычисления обёрнута в счётчик, так что числа в выводе — честные числа полётов.

У троичного поиска есть одно удобство: его трудно написать неправильно. В золотом сечении больше переменных, и легко перепутать, какая точка куда переходит. Но если вычисление функции дорогое, как здесь, лишние 25 полётов окупают любую аккуратность.

Если известен наклон

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

Самый частый случай — функция на целых числах. Вместо производной сравнивают f(m) и f(m + 1): если следующее значение больше, вершина правее m. Это двоичный поиск первого m, где f(m) ≥ f(m + 1), и он короче и надёжнее троичного поиска на целых с его концовкой «осталось три точки — проверить все». Дихотомия из листинга — то же самое для вещественных чисел, где роль «плюс один» играет крошечный сдвиг.

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

Вершину можно найти только наполовину

У поиска вершины есть ограничение, которого нет у поиска корня. Рядом с вершиной функция плоская: f(x* + d) ≈ f(x*) − c·d². Число double хранит около шестнадцати значащих цифр, и как только c·d² становится меньше последней из них, две точки неразличимы. Значит, положение вершины известно лишь с точностью около √ε — восьми знаков из шестнадцати, сколько ни проси.

   asked   golden: |x - 1/2|   bisection of f′
   1e-08             0.0e+00           3.7e-09
   1e-10             6.4e-10           2.9e-11
   1e-12             6.8e-10           4.5e-13
   1e-14             6.8e-10           3.6e-15

Золотое сечение застревает на 6,8·10⁻¹⁰: дальше точки неразличимы, и ответ решает округление. Двоичный поиск по знаку производной 1 − 2x доходит до 3,6·10⁻¹⁵: у корня функция не плоская, а крутая. Поэтому библиотечные функции не принимают точность меньше примерно √ε: у optimize в R точность по умолчанию — корень четвёртой степени из машинного эпсилона, около 1,2·10⁻⁴. Хорошая новость в том же: если вершина плоская, то и значение функции в найденной точке почти не отличается от максимума. Для броска ошибка угла в тысячную градуса стоит меньше миллиметра дальности.

Где ломается: шум, площадки, несколько вершин

Всё сказанное держится на одном предположении: функция унимодальна. У симуляции полёта это так. У настоящих данных — почти никогда. Вот средняя суточная температура в Алматы за каждый год с 1941-го. Летом жарко, зимой холодно: кажется, вершина одна. Попробуйте найти самый жаркий день года троичным поиском.

В 1979 году самым жарким был день 3 августа, 28,3 °C. Троичный поиск называет 24 июня, 22,2 °C. На четвёртом шаге он сравнил два дня в середине лета, второй оказался прохладнее, и вся вторая половина лета, вместе с настоящей жарой, ушла в отбросы. Дальше поиск честно сходится к вершине маленького холма в июне. Из 85 лет он ошибается в 53. Сглаживание помогает, но не до конца: по неделям промахов 36, по месяцам — 13, по двухмесячному окну — всё ещё 10.

Причина не в данных, а в предположении. День в середине лета может оказаться холоднее соседних из-за дождя, и тогда у кривой не одна вершина, а десятки. Сравнение двух точек отвечает на вопрос «где выше в этих двух местах», но не «где вершина». Для поиска максимума в данных есть простой и надёжный способ: пройти их все. 365 сравнений — это меньше миллисекунды.

Отдельная беда — площадки. Если функция где-то постоянна, равенство f(m₁) = f(m₂) ничего не говорит о том, где вершина, и никакая хитрость не поможет.

Так что на целых числах с повторами — скажем, у функции, которая округляет время до минуты, — троичный поиск может работать, а может и нет, в зависимости от того, где площадки. Строгая унимодальность или выпуклость нужна по-настоящему. Выпуклые функции — самый надёжный источник унимодальных: сумма выпуклых выпукла, максимум выпуклых выпуклый. На этом построена вторая задача статьи.

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

Сколько это стоит

Время — число вычислений функции: O(log(L/ε)), где L — длина отрезка, ε — нужная точность. Константа зависит от способа: на одно вычисление троичный поиск сжимает отрезок в 1,22 раза, золотое сечение — в 1,62. Памяти — несколько переменных. На целых числах от 1 до n — O(log n) вычислений.

Замер ниже — время поиска угла для теннисного мяча при разной требуемой точности: от десятой доли градуса до 10⁻¹⁰. Каждый полёт стоит одинаково, поэтому прямые идут с постоянным отношением: троичный поиск медленнее золотого сечения примерно в 2,3 раза при любой точности.

Время растёт линейно по числу знаков: каждый новый знак стоит одинакового числа полётов. Обратите внимание, что за восьмым знаком оно растёт впустую: как показано выше, вершина дальше не уточняется.

Где работает

Библиотеки оптимизации. В SciPy функция minimize_scalar умеет искать золотым сечением, но по умолчанию берёт метод Брента: золотое сечение, которое, где можно, ускоряется параболой через три известные точки. Документация так и советует: Брент обычно лучше. Функция optimize в R — тот же метод Брента по его книге 1973 года «Algorithms for Minimization without Derivatives».

Выбор шага в многомерной оптимизации. Метод, который выбрал направление спуска, должен решить, как далеко идти. Это поиск минимума функции одной переменной — длины шага, и в учебниках его решают золотым сечением или методом Брента.

Олимпиадные задачи. Минимум выпуклой функции: расстояние между движущимися точками, разброс группы, стоимость, которая складывается из растущей и убывающей частей. Если функция зависит от двух переменных и выпукла, поиск вкладывают: для каждого x внутренним поиском находят лучший y. Это законно, потому что минимум выпуклой функции по одной переменной — выпуклая функция другой.

Где не работает. Поиск максимума в массиве данных, как в сцене с температурой: там нужен полный просмотр. И поиск в отсортированном массиве: для него есть двоичный поиск, а троичный, который иногда предлагают «потому что три больше двух», делает больше сравнений: две проверки отбрасывают лишь две трети.

История

В 1953 году статистик Джек Кифер опубликовал заметку «Sequential minimax search for a maximum». Он искал способ найти максимум унимодальной функции за заданное число вычислений так, чтобы в худшем случае остаток отрезка был как можно короче, — и доказал, что лучший способ — поиск Фибоначчи. Золотое сечение получается из него, когда вычислений много. В 1973 году Ричард Брент в книге об оптимизации без производных соединил золотое сечение с параболической интерполяцией: парабола быстро находит вершину гладкой функции, а золотое сечение страхует там, где парабола ошибается.

Задача о броске старше самих компьютеров. Что в пустоте дальше всего летит снаряд, брошенный под 45°, знали со времён Галилея. Сопротивление воздуха, которое опускает этот угол, мучило баллистиков веками, и одной из первых задач электронных вычислительных машин были именно таблицы стрельбы — те самые полёты шаг за шагом, что в нашей сцене. Настоящие толкатели ядра выпускают его ниже теоретического угла по другой причине: чем круче выпуск, тем меньше скорость, которую человек успевает передать ядру. Ник Линторн в 2001 году показал на пяти спортсменах, что у каждого свой оптимум и он совпадает с привычным спортсмену углом.

Задачи

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