Алгоритм Флойда — Уоршелла

Кокшетау и Петропавловск стоят в 163 километрах друг от друга. А кратчайший путь между ними по воздуху — 2628 километров: из Кокшетау в Алматы и оттуда обратно на север. Прямых рейсов между ними в базе нет, и самолёт обходит пол-карты. Таких ответов нужно 420 — по одному на каждую упорядоченную пару из 21 аэропорта, — и вместе они складываются в таблицу. Ниже её заполняет алгоритм Флойда — Уоршелла: на каждом шаге он разрешает пересадку ещё в одном аэропорту и проходит по таблице строка за строкой. Коснитесь клетки, чтобы увидеть маршрут, и аэропорта на карте, чтобы закрыть его.

Сначала в таблице только прямые рейсы — 108 пар из 420. В порядке «сначала узлы» первым шагом разрешается пересадка в Алматы, и маршрут появляется у 382 пар, после Астаны — у всех. Окончательной таблица становится уже после седьмого шага: дальше аэропорты ничего не улучшают. Переключите порядок на «узлы в конце» — и тринадцать шагов подряд не изменится ни одна клетка, а всё решат последние два. Таблица в конце получается одна и та же, а вот путь к ней разный: нижняя полоска показывает, сколько клеток уже окончательны после каждого шага.

Маршруты берутся из открытой базы Virtual Radar Server. Её пополняют пользователи этой программы, которые принимают сигналы самолётов: для каждого позывного там записано, откуда и куда летает рейс. Это не расписание, часть рейсов может быть сезонной или устаревшей, поэтому таблица отвечает на вопрос «как долететь по рейсам этой базы», а длина маршрута — сумма расстояний по дуге большого круга. По земле эти 163 километра, конечно, проезжают без всякого Алматы.

Идея

Кратчайший маршрут из i в j — либо прямой рейс, либо цепочка с пересадками. Перебирать цепочки бессмысленно, их слишком много. Флойд — Уоршелл перебирает другое: не маршруты, а множество аэропортов, где пересаживаться разрешено. Пронумеруем аэропорты и обозначим через d_k[i][j] длину кратчайшего маршрута из i в j, все пересадки которого — в аэропортах с номерами меньше k. При k = 0 пересаживаться нельзя нигде, и d₀ — таблица прямых рейсов. Когда разрешены все, d_n — ответ.

Разрешим пересадку ещё в одном аэропорту k. Лучший маршрут из i в j либо через k не идёт — тогда он уже записан в таблице, — либо идёт через k ровно один раз, и тогда он складывается из лучшего маршрута от i до k и лучшего от k до j, а оба уже в таблице. Значит, новая таблица получается из старой одной формулой для каждой клетки.

Это динамическое программирование в чистом виде: подзадача — «маршрут с пересадками только в первых k аэропортах», переход — один новый аэропорт. В коде k идёт во внешнем цикле, а i и j пробегают всю таблицу:

for k in range(n):
    for i in range(n):
        for j in range(n):
            d[i][j] = min(d[i][j], d[i][k] + d[k][j])

В сцене это видно буквально. Строка и столбец аэропорта k подсвечены: на шаге k каждая строка i берёт из своего столбца k одно число — сколько лететь до k — прибавляет его ко всей строке k и оставляет в каждой клетке меньшее из двух.

Почему таблица верна

Остаётся вопрос, который любят задавать на собеседованиях: формула говорит о двух таблицах, d_k и d_k+1, а код пишет в одну. Не испортит ли новое число в d[i][k] следующие клетки строки? Не испортит, потому что строка k и столбец k на шаге k не меняются: d[i][k] = min(d[i][k], d[i][k] + d[k][k]), а d[k][k] = 0. Все клетки шага читают из строки и столбца k одни и те же числа, в каком порядке их ни обходи. Поэтому памяти нужно O(n²), а не O(n³).

Ловушка: порядок циклов

Три цикла выглядят симметрично, и рука сама пишет их в привычном порядке: for i, for j, for k. Тогда каждая клетка (i, j) считается один раз, перебором всех пересадок k, — кажется, даже естественнее. Но доказательство выше держалось на том, что к моменту шага k клетки d[i][k] и d[k][j] уже окончательны для всех более ранних пересадок. В порядке i, j, k это не так. Вот шесть аэропортов запада, пронумерованных так же, как в данных, — по коду ИКАО. Сцена начинается с конца первого прохода неправильного цикла.

Строка Туркестана считается первой. Когда очередь доходит до клетки «Туркестан → Актау», клетка «Туркестан → Актобе» ещё пуста — её очередь в этой же строке позже, — а строку Костаная ещё не считали вовсе. А настоящий маршрут идёт как раз через Костанай и Актобе: 1160 + 547 + 847 = 2554 км. Проход заканчивается, и три клетки остаются пустыми, хотя маршруты есть. Правильному порядку на ту же таблицу хватает 27 сравнений — считаем только те, где обе половины маршрута уже известны, — а неправильный ошибается после 31 сравнения, и всё исправляет только второй проход.

На всей сети ошибка тише: узлы спасают почти всё, и неправильный цикл ошибается в двух клетках из 420 — из Туркестана в Актау он насчитал 2604 км через Астану вместо 2554 через Костанай и Актобе. Порядок i, k, j ошибается в семи клетках. Именно поэтому ошибку легко пропустить: на маленьком тесте ответ сходится, на большом почти сходится.

Неправильный цикл можно спасти повторением. В 2019 году Икуми Хидэ, Со Кумабэ и Такэнори Маэхара доказали, что порядки i, j, k и i, k, j, повторённые три раза подряд, всегда дают верную таблицу, — статья так и называется: «неправильные реализации Флойда — Уоршелла дают верный ответ после трёх повторов». Двух бывает мало: путь 1 → 5 → 3 → 2 → 6 → 4 → 0 из листинга после второго прохода всё ещё не знает дороги из 1 в 0. Но втрое больше работы ради ошибки в одной строке — плохой обмен. Запомнить проще правило: k — снаружи.

Код

Кроме расстояний, программа хранит вторую таблицу — nxt[i][j], первый перелёт лучшего маршрута из i в j. Когда маршрут через k оказывается короче, его первый перелёт — тот же, что у маршрута из i в k, поэтому nxt[i][j] = nxt[i][k]. Чтобы восстановить весь маршрут, идём по этим ссылкам от i, пока не придём в j. Таблица одна на все пары, и путь между любыми двумя аэропортами читается за число его перелётов.

Прямой рейс есть у 108 пар, одной пересадкой обходятся 299, двумя — 13: дальше трёх перелётов не бывает. Крюки бывают огромными: из Астаны в Караганду, до которой 200 км, кратчайший маршрут по рейсам базы идёт через Алматы и тянется на 1708 км. Самый короткий из самых дальних путей — у Шымкента: из него до любого аэропорта не больше 1935 км. Такую вершину называют центром графа, и это не Алматы и не Астана: у Алматы самый дальний путь — 2097 км, у Кокшетау — 3023 км.

Последние три строки — сцена в цифрах: сколько клеток окончательны после каждого шага при трёх порядках. В Python внутренний цикл ускоряют двумя приёмами: строки d[i] и d[k] берут в локальные переменные, а строку, где d[i][k] бесконечно, пропускают целиком — через k из неё всё равно не долететь. Ещё быстрее отдать весь шаг k numpy, как в главе курса: np.minimum(d, d[:, k, None] + d[None, k, :], out=d).

Тот же цикл без чисел: кто от кого зависит

Заменим числа на «да» и «нет»: reach[i][j] — можно ли попасть из i в j. Минимум превращается в «или», сумма — в «и», и формула шага читается так: из i можно попасть в j, если можно было и раньше или если можно попасть в k, а из k — в j. Это и есть алгоритм Уоршелла 1962 года — он вычисляет транзитивное замыкание. Строку таблицы удобно хранить битами: тогда «если из i есть путь в k, добавить в строку i всю строку k» — одна операция OR над машинными словами, по 64 клетки за раз.

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

Пакету fmt нужно 58 пакетов, net/http — 186, и команда go list -deps печатает ровно столько же плюс сам пакет. От unsafe зависит 331 пакет из 379. Порядок k и здесь меняет только дорогу: если идти от нижних пакетов к верхним, к пятидесятому шагу известно 12 218 пар из 21 504, а в обратном порядке — 2980.

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

Отрицательные рёбра и циклы

Доказательство требовало одного: циклы неотрицательны. Отрицательные рёбра ему не мешают, и в этом отличие от Дейкстры. Если же отрицательный цикл есть, кратчайших путей через него не существует, а алгоритм это замечает сам: на диагонали появляется отрицательное число — путь из i обратно в i короче нуля. Поэтому после алгоритма достаточно проверить d[i][i] < 0. Сам такой цикл и его смысл — круг обменов, который возвращает больше, чем было, — разобраны в статье о Беллмане — Форде. При отрицательном цикле числа в таблице могут расти по модулю очень быстро, и на целых фиксированной длины их стоит обрывать, как только диагональ ушла в минус.

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

Три цикла по n — O(n³) времени, таблица — O(n²) памяти, и число рёбер не влияет ни на то, ни на другое. Альтернатива — Дейкстра из каждой вершины: O(n (n + m) log n) с кучей. Для разреженного графа, где m порядка n, это O(n² log n), заметно лучше. Для плотного, где m порядка n², — O(n³ log n), хуже. Отрицательные рёбра Дейкстра не умеет, но их можно убрать потенциалами — это алгоритм Джонсона. Замер на случайных графах: четыре ребра на вершину (редкий граф) и n рёбер на вершину (плотный).

Наклон Флойда — Уоршелла на Go близок к 3: вдвое больше вершин — в восемь раз больше времени; на Python наклон в замере меньше, около 2,3. На Python в редком графе на 200 вершинах Дейкстра из каждой вершины быстрее почти в 10 раз, а в плотном на 100 вершинах они почти сравнялись: Флойд — Уоршелл медленнее всего в 1,4 раза. В Go граница сдвигается: у Флойда — Уоршелла нет кучи и указателей, только проход по строкам подряд, и даже на редком графе он быстрее до четырёхсот вершин, а в плотном на 400 вершинах — в 6 раз. Версия на numpy на тех же 200 вершинах быстрее чистого Python в 42 раза: шаг k — одна операция над всей таблицей.

Отсюда правило выбора. Нужны все пары, а вершин сотни или пара тысяч — Флойд — Уоршелл: он короче, не ошибается на отрицательных рёбрах и хорошо ложится на векторные операции и кэш. Вершин десятки тысяч и больше — таблица n² сама по себе не поместится в память, и нужны пары по запросу, то есть Дейкстра, A* или предподсчёт вроде иерархий сжатия. Для улиц Алматы из 7909 перекрёстков таблица заняла бы 62,5 миллиона клеток.

Где работает

• Библиотеки. В SciPy есть scipy.sparse.csgraph.floyd_warshall, в NetworkX — floyd_warshall и floyd_warshall_numpy: когда нужна вся таблица расстояний небольшого графа, берут их.
• Замыкания отношений. Кто от кого зависит, кто кому подчиняется, какие состояния достижимы — везде, где отношение задано матрицей и нужно его замыкание. Битовая версия обрабатывает 64 пары за операцию.
• Автоматы и регулярные выражения. Алгоритм Клини, по которому из автомата строят регулярное выражение, — тот же тройной цикл по «разрешённым промежуточным состояниям».
• Системы реального времени. В проверке моделей с таймерами ограничения вида «x − y ≤ c» хранят матрицей, и чтобы привести её к каноническому виду с самыми тесными ограничениями, по ней проходят Флойдом — Уоршеллом.

Три статьи 1959–1962 годов

В 1959 году Бернар Руа напечатал в докладах Академии наук заметку «Транзитивность и связность» — с тем же способом вычислять замыкание. В январе 1962 года в Journal of the ACM вышла «Теорема о булевых матрицах» Стивена Уоршелла. Как рассказывают, Уоршелл поспорил с коллегой на бутылку рома, кто первым выяснит, всегда ли этот алгоритм работает, — и принёс доказательство наутро; ромом, говорят, поделился с проигравшим. В июне того же года в Communications of the ACM Роберт Флойд опубликовал «Algorithm 97: Shortest Path» — одну страницу с той же схемой для кратчайших путей. В том же году Питер Ингерман в «Algorithm 141: Path Matrix» записал замыкание тремя вложенными циклами — в том виде, в каком его пишут сейчас.

Отсюда и разнобой в названиях: алгоритм Флойда, Уоршелла, Роя — Уоршелла, Роя — Флойда. А сама идея — перебирать не пути, а множество разрешённых промежуточных вершин — восходит к работе Стивена Клини 1956 года о событиях в нервных сетях и конечных автоматах, где появились регулярные выражения.

Задачи

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