Просмотр кода страницы Алгоритм поиска кратчайших путей между всеми парами при помощи матричного произведения

Ключевые слова и синонимы: кратчайший путь; анализ алгоритмов

== Постановка задачи ==
Задача поиска кратчайших путей между всеми парами (APSP) заключается в нахождении кратчайших путей между всеми парами вершин в ориентированном графе с ребрами неотрицательной действительной стоимости. Основное внимание уделяется кратчайшим расстояниям между вершинами, поскольку кратчайшие пути могут быть найдены с незначительным возрастанием стоимости. Классические алгоритмы решают задачу нахождения кратчайших путей между всеми парами вершин за кубическое время – <math>O(n^3) \, </math>. Наша задача заключается в достижении меньшего времени выполнения алгоритма для графа с небольшими целочисленными стоимостями ребер.

Пусть дан ориентированный граф G = (V, E), где V = {1, ..., n} – множество вершин, а E – множество ребер. Стоимость ребра <math>(i, j) \in E\, </math> обозначается как <math>d_{ij}\, </math>. Матрица D размера <math>n \times n</math> содержит элемент <math>d_{ij}\, </math> в позиции (i, j). Для простоты предположим, что <math>d_{ij} > 0 \, </math> и <math>d_{ii} = 0 \, </math> для всех <math>i \ne j</math>. Если не существует ребра из i в j, положим <math>d_{ij} = \infty</math>. Стоимость, или расстояние, пути представляет собой сумму стоимостей всех ребер на пути. Длина пути равна количеству ребер на пути. Кратчайшее расстояние между вершинами i и j равняется минимальной стоимости среди всех путей, ведущих от i до j, и обозначается как <math>d_{ij}^*</math>. Пусть <math>D^* = {d_{ij}^*}</math>. Значение n называется размером матрицы.

Пусть A и B – матрицы (n, n). Умножение матриц A и B может производиться одним из трех способов:
(1) Обычное матричное произведение над кольцом <math>C = AB \, </math>; (2) булево произведение матриц <math>C = A \cdot B</math>; (3) произведение матриц расстояния <math>С = A \times B</math>, где

(1) <math>c_{ij} = \sum_{k=1}^n a_{ik} b_{kj}</math>

(2) <math>c_{ij} = \bigvee_{k=1}^n a_{ik} \wedge b_{kj}</math>

(3) [[Файл:Cij_min.jpg]]

В любом из трех случаев матрица C называется произведением, а процесс вычисления – умножением. Во всех трех случаях k пробегает весь диапазон (1, ..., n}. Частичное произведение матриц A и В определяется для значений k из подмножества I множества V. Иными словами, частичное произведение получается при умножении вертикальной прямоугольной матрицы A(*, I), столбцы которой извлечены из матрицы A согласно множеству I, и горизонтальной прямоугольной матрицы B(I, *), полученной из B путем извлечения строк, соответствующих множеству I. Интуитивно множество I представляет собой множество контрольных точек k при переходе от i к j в графе.

Лучший алгоритм [3] вычисляет произведение (1) за время <math>O(n^ \omega) \, </math>, где <math>\omega \, </math> = 2.376. Далее мы будем учитывать три знака после запятой. Для вычисления произведения (2) булевы значения 0 и 1 в матрицах A и B можно рассматривать как целые числа, применить алгоритм для вычисления (1) и затем заменить ненулевые элементы получившейся матрицы на 1. Следовательно, его сложность также равна <math>O(n^ \omega) \, </math>. ''Свидетели'' произведения (2) задаются матрицей свидетелей <math>W = {w_{ij}} \, </math>, где <math>w_{ij} = k \, </math> для некоторого k, такого, что <math>a_{ik} \wedge b_{kj} = 1</math>. Если такого k не существует, <math>w_{ij} = 0 \, </math>. ''Матрица свидетелей'' <math>W = w_{ij} \, </math> для (3) определяется следующим образом: <math>w_{ij} = k \, </math>, при котором значение <math>c_{ij} \,</math>  минимально. Если существует алгоритм для (3) с временем T(n), игнорирующий полилогарифмический коэффициент n, задача поиска кратчайших путей между всеми парами может быть решена за время Õ(T(n)) путем последовательного применения метода возведения в квадрат, заключающегося в повторном выполнении <math>D \leftarrow D \times D</math> <math>O(\log n) \,</math> раз.

Вычисление кратчайших путей будет здесь определяться как нахождение матрицы свидетелей размера n, которой может задаваться кратчайший путь из i в j, за время O(ℓ), где ℓ – длина пути. Говоря более точно, если <math>w_{ij} = k \, </math> в матрице свидетелей <math>W = w_{ij} \, </math>, это означает, что путь из i в j проходит через k. Таким образом, рекурсивная функция path(i,j) задается равной <math>(path(i,k),k, path(k, j))\, </math>, если <math>w_{ij} = k > 0 \, </math>, и <math>nil\, </math> в случае <math>w_{ij} = 0 \, </math>, когда путь определяется списком вершин, исключая конечные точки. В дальнейшем изложении k будет записываться в <math>w_{ij} \, </math> в тех случаях, когда будет найден такой элемент k, что путь из i в j модифицируется или задается заново путями из i в k и из k в j. Предшествующие результаты приводятся в качестве схемы получения основных результатов.

== Алгоритм Алона, Галила и Маргалита ==

Полное изложение алгоритма приведено в [1]. Пусть все ребра графа имеют единичную стоимость. Пусть <math>D^{(\ell)} </math> – <math>\ell </math>-я приближенная матрица для <math>D^* \, </math>, определенная следующим образом: <math>d_{ij}^{(\ell)} = d^*_{ij}</math>, если <math>d^*_{ij} ≤ \ell \, </math>, и <math>d_{ij}^{(\ell)} = \infty \, </math> в противном случае. Пусть A – матрица смежности графа G: <math>a_{ij} = 1 \, </math>, если существует ребро (i, j), и <math>a_{ij} = 0 \, </math> в противном случае. Пусть <math>a_{ii} = 1 \, </math> для всех i. Алгоритм состоит из двух этапов. На первом этапе вычисляются <math>D^{(\ell)} </math> для <math>\ell = 1, ..., r</math> при помощи проверки (i, j)-го элемента <math>A^{\ell} = {a^{\ell}_{ij}}</math>. Заметим, что если <math>a^{\ell} = 1</math>, то существует путь из i в j длиной <math>\ell</math> или меньше. Поскольку булево произведение матриц может быть вычислено за время <math>O(n^{\omega}\, )</math>, время выполнения этого этапа составит <math>O(rn^{\omega}\, )</math>.

На втором этапе алгоритм вычисляет <math>D^{(\ell)} </math> для <math>\ell = r, \mathcal{d} 3/2 r \mathcal{e}, \mathcal{d} 3/2 \mathcal{d} 3/2 r \mathcal{e}\mathcal{e}, ... , n’</math> путем последовательного возведения в квадрат, где n’ – наименьшее целое число в последовательности <math>\ell \,</math>, такое, что <math>\ell \ge n</math>. Пусть <math>T_{i\alpha} = \, </math> {<math>j \mid d_{ij}^{(\ell)} = \alpha</math>} и <math>I_i = T_{i\alpha} \, </math>, такое, что <math>\mid T_{j\alpha}\mid</math> минимально для <math>\mathcal{d} \ell /2 \mathcal{e} < \alpha < \ell</math>. Основное наблюдение второго этапа заключается в следующем: необходимо проверять k только в I<math>_i</math>, размер которого не превышает <math>2n/ \ell \, </math>, поскольку корректные расстояния между <math>\ell + 1 \, </math> и <math>\mathcal{d} 3 \ell/2 \mathcal{e}</math> могут быть получены как суммы <math>d_{ik}^{(\ell)} + d_{kj}^{(\ell)}</math> для некоторых k, удовлетворяющих неравенству <math>\mathcal{d} \ell /2 \mathcal{e} \le d_{ik}^{(\ell)} \le \ell</math>. Значение I<math>_i</math> сходно с I для частичных произведений – за тем исключением, что I меняется для каждого i. Таким образом, время одного возведения в квадрат составляет <math>O(n^3/ \ell \,)</math>; а полное время второго этапа задается при <math>N = \mathcal{d}log_{\frac{3}{2}} n/r\mathcal{e}</math> формулой <math>O(\sum_{s=1}^N n^3/((3/2)^s r) = O(n^3/r)</math>. Сбалансировав оба этапа условием <math>rn^{\omega} = n^3/r \,</math>, получаем время <math>O(n^{(\omega + 3)/2}) = O(n^{2.688}) \,</math> для алгоритма с <math>r = O(n^{(3- \omega)/2}) \, </math>.

Свидетели могут быть вычислены на первом этапе за время, полилогарифмическое относительно n, по методу, приведенному в [2]. Поддержка свидетелей на втором этапе тривиальна.
Если дан ориентированный граф G, стоимости ребер которого задаются целыми числами от 1 до M, где M – целое положительное число, то граф G можно преобразовать в G’, заменив каждое ребро соответствующим количеством ребер единичной стоимости, вплоть до M. Очевидно, что задача для графа G может быть решена путем применения вышеприведенного алгоритма к G’, для чего потребуется время <math>O((Mn)^{(\omega +3)/2} \, )</math>. Это время оказывается меньше кубического при <math>M < n^{0.116} \, </math>. Поддержка свидетелей вносит дополнительный полилогарифмический коэффициент в каждом случае.

Для неориентированных графов с единичной стоимостью ребра Зейделем было доказано время выполнения <math>\tilde{O}(n^{\omega}) \, </math> [7].

== Алгоритм Такаоки ==

Если стоимости ребер ограничены положительным целым числом M, можно разработать более удачный алгоритм, что было показано Такаокой в [9]. Вкратце рассмотрим алгоритм Романи [6] для произведения матриц расстояния.

Пусть A и B – матрицы расстояния формата (n, m) и (m, n), соответственно, элементы которых ограничены числом M либо равны бесконечности. Пусть диагональные элементы равны 0. Преобразуем матрицы A и B в A' и B', где <math>a'_{ij} = (m + 1)^{M - a_{ij}}</math>, если <math>a_{ij} \ne \infty \, </math>, и 0 в противном случае, а <math>b'_{ij} =(m + 1)^{M - b_{ij}}</math>, если <math>b_{ij} \ne \infty \, </math>, и 0 в ином случае.

Пусть C' = A' B' – результат обычного матричного произведения, а C = A × B – результат произведения матриц расстояния. Тогда верно следующее:

<math>c'_{ij} = \sum_{k=1}^m (m + 1)^{2M - (a_{ik} + b_{kj})}</math>, <math>c_{ij} = 2M - \mathcal{b} log_{m + 1} c'_{ij} \mathcal{c}</math>

Такое произведение матриц расстояния называется произведением (n, m)-Романи. Ранее мы имели дело с умножением квадратных матриц, или произведением (n, n)-Романи. В следующем разделе будет рассмотрен случай для m < n.

Матрица С может быть вычислена при помощи <math>O(n^{\omega}) \, </math> арифметических операций на целых числах вплоть до <math>(n + 1)^M \, </math>. Поскольку эти значения можно выразить при помощи O(M log n) бит, а алгоритм Шонхаге и Штрассена [8] для умножения k-битных чисел требует O(k log k log log k) битовых операций, для вычисления C потребуется время O(<math>n^{\omega}</math> M log n log(M log n) log log(M log n)), или <math>\tilde{O}(Mn^{\omega}) \, </math>.

Первый этап заменяется алгоритмом на основе произведения (n, n)-Романи, а второй этап модифицируется таким образом, чтобы иметь дело с длинами путей, а не расстояниями.

Заметим, что граница M была заменена на <math>\ell</math>M в произведении матриц расстояния на первом этапе. Игнорируя полилогарифмические коэффициенты, получаем время выполнения для первого этапа – <math>\tilde{O}(n^{\omega} r^2 M) \, </math>. Мы предполагаем, что M соответствует <math>O(n^k) \, </math> для некоторого константного k. Сравнивая эту сложность со сложностью второго этапа, <math>O(n^3/r) \, </math>, получаем полное время вычисления <math>\tilde{O}(n^{(6 + \omega)/3} M^{1/3}) \, </math> при выборе <math>r = O(n^{(3 - \omega)/3} M^{-1/3}) \, </math>. Чтобы время выполнения оказывалось меньше кубического, значение M должно почти совпадать с <math>O(n^{0.624}) \, </math>.

== Основные результаты ==

Цвик улучшил алгоритм Алона-Галила-Маргалита в нескольких отношениях. Самым серьезным улучшением стало снижение времени нахождения кратчайших путей между всеми парами с <math>O(n^{2.688}) \, </math> до <math>O(n^{2.575}) \, </math> для ребер единичной стоимости. Наибольшее ускорение в алгоритме Алона-Галила-Маргалита [1] было достигнуто в быстрой операции булева произведения матриц, а в алгоритме Такаоки [9] – в быстром умножении матриц расстояния по Романи; в обоих случаях мы имели дело с быстрым матричным умножением квадратных матриц.
 
В данном разделе «ускорителем» служит алгоритм быстрого умножения матриц расстояния по Романи, созданный на основе алгоритма быстрого умножения прямоугольных матриц, разработанного Копперсмитом [4] и Хуангом и Пэном [5]. Пусть <math>\omega(p, q, r) \, </math> – порядок временной сложности при умножении матриц <math>(n^p, n^q) \, </math> и <math>(n^q, n^r) \, </math>. Предположим, что произведение матриц (n, m) и (m, n) может быть вычислено за <math>O(n^{\omega (1, \mu, 1)}) \, </math> арифметических операций, где <math>m = n^{\mu} \, </math> и <math>0 \le \mu \le 1 \, </math>. Нам известно следующее: <math>O(n^{\omega( 1, 1, 1)}) = O(n^{2.376}) \, </math> и <math>O(n^{\omega (1, 0.294, 1)}) = \tilde{O}(n^2) \, </math>. Для вычисления произведения квадратных матриц (n, n) требуется <math>n^{1 - \mu} \, </math> операций матричного умножения, что дает время выполнения <math>O(n^{(\omega(1, \mu, 1) + 1 - \mu}) \, </math>, которое можно преобразовать в <math>O(n^{2 + \mu}) \, </math>, где <math>\mu \, </math> удовлетворяет равенству <math>\omega(1, \mu, 1) = 2 \mu + 1 \, </math>. На данный момент лучшее известное значение <math>\mu \, </math> составляет <math>\mu = 0.575 \,</math>, так что временная граница превращается в <math>O(n^{2.575}) \, </math> – это значение уступает <math>O(n^{2.376}) \, </math>. Поэтому в дальнейшем будет использоваться алгоритм для прямоугольных матриц.

Вышеприведенный алгоритм встроен в вычисление произведения (n, m)-Романи при <math>m = n^{\mu} \, </math> и <math>M = n^t \, </math> для некоторого t > 0, а время выполнения составляет <math>\tilde{O}(Mn^{\omega (1, \mu, 1)}) \, </math>. Следующий этап заключается во встраивании произведения (n, m)-Романи в алгоритм нахождения кратчайших путей между всеми парами. Первый алгоритм представляет собой последовательное применение операции возведения в квадрат, напоминающее второй этап алгоритма в [1]. Чтобы с пользой применить алгоритм произведения прямоугольных матриц (n, m)-Романи, нам потребуется определение мостового множества, которое будет играть ту же роль, что и множество I в алгоритме частичного вычисления матричного произведения в разделе «Постановка задачи».

Обозначим за <math>\delta(i,j) \, </math> кратчайшее расстояние между i и j, а за <math>n(i,j) \, </math> – минимальную длину среди всех кратчайших путей от i до j. Подмножество I множества вершин V является <math>\ell</math>-мостовым множеством, если удовлетворяет следующему условию: если <math>\eta(i,j) \ge \ell</math>, существует <math>k \in I</math>, такое, что <math>\delta(i,j) = \delta(i,k) + \delta(k,j) \, </math>. I является сильным <math>\ell</math>-мостовым множеством при выполнении условия: если <math>\eta(i,j) \ge \ell \, </math>, то существует <math>k \in I</math>, такое, что <math>\delta(i,j) = \delta(i,k) + \delta(k,j) \, </math> и <math>\eta(i,j) = \eta(i,k) + \eta(k,j) \, </math>. Заметим, что для графа с ребрами единичной стоимости эти множества совпадают.

Заметим также, что если <math>(2/3) \ell \le \mu(i,j) \le \ell \,</math> и I является сильным <math>\ell/3</math>-мостовым множеством, существует <math>k \in I</math>, такое, что <math>\delta(i,j) = \delta(i,k) + \delta(k,j) \, </math> и <math>\mu(i,j) = \mu(i, k) + \mu(k,j) \,</math>. Если выполняется свойство сильного мостового множества, произведение (n, m)-Романи может быть использовано для решения задачи нахождения кратчайших путей между всеми парами следующим образом. Последовательным применением возведения в квадрат, как в случае алгоритма Алона-Галила-Маргалита, алгоритм вычисляет <math>D^{(\ell)} \, </math> для <math>\ell = 1, \mathcal{d} 3/2 \mathcal{e}, \mathcal{d} 3/2 \mathcal{d} 3/2 \mathcal{e} \mathcal{e}, ... , n' \, </math> (где n' – первое значение <math>\ell</math>, превышающее n), используя различные варианты множества I, описанного ниже. Для вычисления мостового множества алгоритм поддерживает матрицу свидетелей с дополнительным полилогарифмическим коэффициентом в оценке сложности. В источнике [10] предложены три способа выбора множества I. Пусть <math>|I| = n^r \, </math> для некоторого r, <math>0 \le r \le 1</math>.

(1) Выберем 9n ln n/ℓ вершин из V случайным образом. В этом случае можно показать, что алгоритм решает задачу нахождения кратчайших путей между всеми парами с высокой вероятностью, т.е. с вероятностью 1 - <math>1/n^c</math> для некоторой константы c > 0; можно показать, что она равна 3. Иными словами, I является сильным <math>\ell/3</math>—мостовым множеством с высокой вероятностью. Время T в основном определяется произведением (n,m)-Романи. <math>T = \tilde{O}( \ell Mn^{(\omega (1, r, 1)})</math>, поскольку элементы матрицы могут иметь величину вплоть до <math>\ell M</math>. Из того, что m = O(n ln n/ℓ) = <math>n^r</math>, следует, что <math>\ell = \tilde{O}(n^{1 - r})</math> и, в силу этого, <math>T = O(Mn^{l - r} n^{\omega (1, r, 1)}) \, </math>. В случае M = 1 граница значения r равна <math>\mu</math> = 0.575, и, следовательно, <math>T = O(n^{2.575}) \, </math>. Если <math>M = n^t \ge 1 \, </math>, время оценивается как <math>О(n^{2 + \mu (t)}) \, </math>, где <math>t \le 3 - \omega = 0.624 \, </math>, а <math>\mu = \mu (t) \, </math> удовлетворяет равенству <math>\omega (1, \mu , 1) = 1 + 2 \mu - t \, </math>. Это определено из наиболее известного <math>\omega(1, \mu, 1) \,</math> и значения t. Поскольку результат является корректным с высокой вероятностью, этот алгоритм является рандомизированным.

(2) Рассмотрим случай с единичными стоимостями ребер. В варианте (1) вычисление свидетелей производится дополнительно, то есть не является необходимым, если нужны только кратчайшие расстояния. Чтобы получить такую же сложность для точного, а не рандомизированного алгоритма, вычисление свидетелей становится обязательным. Как было отмечено ранее, поддержка свидетелей (то есть матрицы W) вносит дополнительный полилогарифмический коэффициент, что означает, что анализ может быть сосредоточен на произведении Романи в <math>\tilde{O} </math>-нотации. Более конкретно, I выбирается в качестве <math>\ell / 3</math>-мостового множества, являющегося сильным при единичной стоимости ребер. Чтобы вычислить I как <math>O(\ell)</math>-мостовое множество, необходимо получить вершины кратчайшего пути из i в j для каждого i и j, используя матрицу свидетелей W, за время <math>O(\ell)</math>. После получения этих <math>n^2 \,</math> множеств за время <math>O(\ell n^2)</math> можно вычислить <math>O(\ell)</math>-мостовое множество размера O(n ln n/ℓ) с той же временной сложностью, как показано в [10]. Процесс вычисления мостового множества должен остановиться в <math>\ell = n^{1/2} \, </math>, поскольку после этой границы процесс становится чрезмерно дорогим; в силу этого после данной границы используется одно и то же мостовое множество. Время вычисления до этой границы остается тем же, что и в варианте (1), а после границы составляет <math>\tilde{O}(n^{2.5}) \, </math>. Таким образом, мы получаем алгоритм из двух этапов.

(3) Если стоимости ребер положительны и ограничены <math>M = n^t > 0 \, </math>, сходная процедура может быть использована для вычисления <math>O(\ell)</math>-мостового множества размера O(n ln n/ℓ) за время <math>\tilde{O}(\ell n^2)</math>. Используя мостовое множество, задачу нахождения кратчайших путей между всеми парами вершин можно решить за время <math>\tilde{O}(n^{2 + \mu (t)}) \, </math> способом, сходным с вариантом (1). Этот результат можно обобщить на случай, когда стоимости ребер находятся в диапазоне от -M до M, сохранив ту же временную сложность, если модифицировать процедуру для вычисления <math>\ell</math>-мостового множества в случае, если не имеется отрицательных контуров. Подробности изложены в [10].

== Применение ==
Эксцентриситетом вершины v называется наибольшее расстояние от v до любой другой вершины графа. Диаметр графа равен наибольшему эксцентриситету для любой вершины. Иными словами, диаметр представляет собой наибольшее расстояние между любой парой вершин. Если задача нахождения кратчайших путей между всеми парами вершин решена, максимальный элемент матрицы будет равен ее диаметру.

== Открытые вопросы ==
Среди оставшихся нерешенными задач особенно выделяются две. Первая заключается в уменьшении сложности <math>\tilde{O}(n^{2.575}) \, </math> при нахождении кратчайших путей между всеми парами вершин графа с ребрами единичной стоимости. Вторая состоит в улучшении границы <math>M < O(n^{0.624}) \, </math> сложности алгоритма нахождения кратчайших путей между всеми парами вершин графа, стоимости ребер которого являются целыми числами не более M, и достижении значений ниже кубических.

== См. также ==

* ''[[Алгоритм поиска кратчайших путей между всеми парами в разреженных графах]]'',

* ''[[Полностью динамический алгоритм нахождения кратчайших путей между всеми парами]]''

== Литература ==

1. Alon, N., Galil, Z., Margalit, O.: On the exponent of the all pairs shortest path problem. In: Proc. 32th IEEE FOCS, pp. 569-575.
IEEE Computer Society, Los Alamitos, USA (1991). Also JCSS 54, 255-262(1997)

2. Alon, N., Galil, Z., Margalit, O., Naor, M.: Witnesses for Boolean matrix multiplication and for shortest paths. In: Proc. 33th IEEE FOCS, pp. 417-426. IEEE Computer Society, Los Alamitos, USA (1992)

3. Coppersmith, D., Winograd, S.: Matrix multiplication via arithmetic progressions. J. Symb. Comput. 9, 251-280 (1990)

4. Coppersmith, D.: Rectangular matrix multiplication revisited. J. Complex. 13,42-49(1997)

5. Huang, X., Pan, V.Y.: Fast rectangular matrix multiplications and applications. J. Complex. 14, 257-299 (1998)

6. Romani, F.: Shortest-path problem is not harder than matrix multiplications. Info. Proc. Lett. 11,134-136 (1980)

7. Seidel, R.: On the all-pairs-shortest-path problem. In: Proc. 24th ACM STOC pp. 745-749. Association for Computing Machinery, New York, USA (1992) Also JCSS 51,400-403 (1995)

8. Schonhage, A., Strassen, V.: Schnelle Multiplikation GroRerZahlen. Computing 7, 281-292 (1971)

9. Takaoka, T.: Sub-cubic time algorithms for the all pairs shortest path problem. Algorithmica 20, 309-318 (1998)

10. Zwick, U.: All pairs shortest paths using bridging sets and rectangular matrix multiplication. J. ACM 49(3), 289-317 (2002)