Ближайшая пара точек

В каталоге USGS 19 073 землетрясения силой от 5 за 2015–2025 годы. Какие два эпицентра ближе всего друг к другу? Сравнить все пары — это 181 880 128 расстояний. Алгоритм ниже обходится тринадцатью с половиной тысячами. Выберите набор точек — землетрясения, аэропорты или города мира — и способ поиска; найденная пара видна на карте и под лупой.

Первый ответ обескураживает: 0 километров. Два землетрясения — 22 апреля 2024 года силой 6,1 и 2 июня 2024 года силой 5,0 — записаны в одной и той же точке с точностью до тысячной доли градуса. Таких совпадений в каталоге три. Одно из них — не два землетрясения, а одно, записанное дважды: в базе USGS у события 27 февраля 2021 года в 18:59:25 силой 5,3 два разных идентификатора при одинаковых времени, месте, глубине и силе. У аэропортов ближайшая пара — запись, которая так и называется: «(Duplicate)YEG», и настоящий аэропорт Эдмонтона. У городов — город и его же район, которые справочник поставил в одну точку.

Так выглядит ближайшая пара почти в любом настоящем каталоге: это дубль. Поэтому её и ищут — при проверке данных первым делом смотрят, какие записи совпадают или почти совпадают. Включите «только разные точки», и ответ станет содержательным: два землетрясения 24 мая 2017 года в 68 метрах друг от друга, второе, сильнее первого, через 29 минут; два аэропорта Красноярска в 3,7 километра.

Идея

Отсортируем точки по x и разрежем вертикальной прямой на две равные половины. Ближайшую пару в каждой половине найдём тем же способом; пусть меньшее из двух расстояний — d. Осталось проверить пары, которые пересекают разрез, а такая пара может быть ближе d, только если обе точки лежат в полосе шириной 2d вокруг разреза. Перебираем точки полосы снизу вверх, и каждую сравниваем только со следующими, пока разница по высоте меньше d: таких точек не больше семи.

Это классический алгоритм Майкла Шеймоса и Дэна Хоуи. Две детали делают его быстрым. Первая — полоса: расстояние до разреза отсекает всё, что дальше d, не считая ни одного расстояния. Вторая — сортировка полосы по y: если сортировать её заново в каждом вызове, работа склейки станет n log n, и всё вместе — n log² n. Вместо этого рекурсия возвращает точки своего куска уже упорядоченными по y, и склейка сливает два упорядоченных списка, как сортировка слиянием, за линейное время.

По шагам

Ниже 174 города от 10 тысяч жителей вокруг Алматы — между 40° и 46° северной широты и 68° и 81° восточной долготы, — положенные на плоскую карту. Запустите поиск: рекурсия режет точки пополам, пока в куске не останется две-три, решает куски и на обратном пути проверяет полосу у каждого разреза. Тонкие линии — посчитанные расстояния, толстая — лучшая пара на этот момент.

Ближайшими оказываются Нижняя Ала-Арча и Маевка, в 1,45 километра друг от друга, и поиск посчитал для этого 168 расстояний из 15 051 возможного. Посмотрите на полосы: пока d большое, в полосу попадает много городов, но по мере того как d уменьшается, полосы сужаются до горстки точек. Переключите «в полосе — все пары», и видно, зачем нужна вторая деталь: без остановки по высоте полосы у верхних разрезов сравниваются каждая с каждой.

Почему в полосе мало кандидатов

Код не считает до семи: он идёт по полосе вверх, пока разница по высоте меньше d, и останавливается. Лемма гарантирует, что таких шагов от каждой точки немного, а проверка по высоте ещё и отсекает точки раньше. Оценка с запасом: подсчитав тщательнее, семь можно уменьшить, но на асимптотику это не влияет.

Код

Классический алгоритм на плоскости. Координаты — целые метры, поэтому сравниваются квадраты расстояний без корней и без ошибок округления. Рекурсия получает номера точек по возрастанию x и возвращает их по возрастанию y; ничьи по координате разрешаются номером точки, чтобы обе версии программы шли одним путём. Кроме городов — два входа, которые выглядят опасно: все точки на одной вертикали, где в полосу попадают все, и решётка, где ничьи на каждом шагу.

one vertical line: 20000 points, closest 6, 15429 distances computed (0.77 per point)
100 x 100 grid: 10000 points, closest 10, 7712 distances computed (0.77 per point)

На вертикали все точки в полосе, но полоса проверяется только до первого соседа выше d, и расстояний выходит меньше одного на точку.

Первая сцена — это та же рекурсия в трёх измерениях. Точка на шаре записана единичным вектором (cos φ cos λ, cos φ sin λ, sin φ), а прямое расстояние сквозь Землю — хорда — растёт вместе с расстоянием по поверхности, так что ближайшие по хорде и есть ближайшие на карте. Разрез по x даёт уже не полосу, а слой; точки слоя перебираются по y, и расстояние считается, только если близка ещё и координата z. Лемма переносится: в слое у каждой точки ограниченное число кандидатов, только константа больше.

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

Сортировка по x один раз — O(n log n). Рекурсия: две половины, склейка — слияние по y, отбор полосы и проход по ней, всё за O(n): T(n) = 2T(n/2) + O(n), второй случай основной теоремы, O(n log n). Память — O(n) на списки и log n кадров рекурсии.

Быстрее в худшем случае нельзя, если алгоритм только сравнивает и складывает координаты: к ближайшей паре сводится вопрос «все ли числа различны» — положите числа на прямую, и они различны, только если ближайшая пара на ненулевом расстоянии. А для него в модели алгебраических деревьев решений известна оценка Ω(n log n).

Замер на случайных точках, разделяй и властвуй против перебора всех пар и против сетки, о которой ниже:

Перебор растёт как n²: на Python уже две тысячи точек — это больше секунды, а 32 тысячи точек заняли бы примерно в 250 раз больше, около пяти минут. Разделяй и властвуй на тех же 32 тысячах — около шести десятых секунды.

Сетка: в среднем O(n)

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

А d заранее неизвестно. Выход нашёл Майкл Рабин в 1976 году, а вариант в листинге — со случайным порядком точек и перестройкой сетки — разбирают учебники, например «Algorithm Design» Джона Клейнберга и Евы Тардош. Берём точки в случайном порядке и держим сетку под текущее лучшее d. Если новая точка дала пару ближе — сетку строим заново под новое d. Перестройка стоит столько, сколько точек уже вставлено, и кажется, что это может обойтись в n². Но в среднем — нет.

Листинг подтверждает оценку: на ста тысячах случайных точек сетка перестраивалась 20 раз и переложила в сумме 185 900 точек — меньше двух на точку. Но «линейно в среднем» не значит «быстрее»: на Python сетка обгоняет рекурсию примерно вдвое, а на Go, где хеш-таблица с ключами-парами обходится дороже рекурсии по срезам, разделяй и властвуй быстрее сетки. Сетка зато проще обобщается: та же раскладка по квадратам отвечает и на вопрос «какие пары ближе r», о котором вторая задача.

Где работает

• Проверка данных. Совпадающие и почти совпадающие координаты — первый признак дублей: одно землетрясение под двумя номерами, один аэропорт двумя записями, город и его район в одной точке. Ближайшая пара отвечает «есть ли дубли вообще», раскладка по сетке — «какие все».
• Контроль сближения. Не подошли ли движущиеся объекты ближе допустимого: в каждый момент это вопрос о ближайшей паре, и так его и решают, когда объектов много.
• Кластеризация. Иерархическая кластеризация одиночной связью начинает с того, что объединяет две ближайшие точки, потом две ближайшие группы.
• Много вопросов к одним точкам. Если нужна не одна пара, а ближайший сосед для многих точек или запросов, строят k-d дерево или триангуляцию Делоне: ближайшая пара всегда соединена ребром триангуляции.

История

Задача о ближайшей паре — одна из первых, с которых в 1970-х началась вычислительная геометрия. Майкл Шеймос и Дэн Хоуи в 1975 году описали алгоритм за O(n log n) на плоскости. В 1976 году Джон Бентли и Шеймос обобщили его на любую фиксированную размерность — с тем же O(n log n), — и в том же году Майкл Рабин показал, что со случайностью и сеткой можно в среднем за O(n); его работа стала одним из первых примеров того, как случайность ускоряет алгоритмы. В 1979 году Стивен Форчун и Джон Хопкрофт получили O(n log log n) без случайности, но с операцией округления вниз, которой нет в модели, где действует нижняя оценка. В 1995 году Самир Хуллер и Йосси Матиас предложили совсем простой вероятностный алгоритм: решето, которое сначала грубо оценивает расстояние, а потом уточняет его сеткой.

Задачи

В первой задаче нужна сама ближайшая пара на настоящих городах мира — и тест, на котором проваливается популярное «сравнить каждую точку с несколькими соседями по x». Во второй — сосчитать все пары ближе заданного расстояния: тут удобнее сетка.