shortest-paths.md

Кратчайшие пути в орграфе

Рассмотрим задачу нахождения кратчайших путей в орграфе от одной вершины ко всем остальным. Мы уже рассматривали похожую задачу в теме BFS, однако там мы считали расстояния в рёбрах. Сейчас рассмотрим так называемые взвешенные графы, то есть графы, каждому ребру которых поставлено в соответствие некоторое число -- вес (или стоимость). Вес можно также воспринимать как функцию c, отображающую множество ребер во множество чисел. Теперь длиной пути будем называть не число ребер, а сумму весов ребер, из которых состоит этот путь. Чтобы облегчить себе задачу, будем предполагать, что веса ребер не могут быть отрицательными.

В связи с наличием дополнительной информации о ребрах изменим наше представление графа. Теперь список смежности вершины v в g[v] будет хранить не список смежных вершин, а список пар смежных вершин и стоимостей (u, c(v, u)).

Будем рассматривать только ориентированные графы, помня, что неориентированные графы обычно представляются в памяти компьютера, как ориентированные, но с продублированными дугами для двух концевых вершин.

Наивный алгоритм

Как и раньше, ответ будем хранить в массиве dist, изначально заполнив его INF. Будем считать, что INF больше любого другого числового значения. Также будем считать, что вершины, которые недостижимы из стартовой вершины start, находятся от нее на растоянии INF.

Считая, что каждая вершина находится на расстоянии 0 от самой себя (c(v, v) = 0), можем сразу же заполнить ответ для стартовой вершины.

dist[start] = 0

Что будем делать дальше? Рассмотрим окружение вершины start. Мы точно знаем, что каждая из вершин u_i окружения N(start) находится не более чем на растоянии c(start, u_i) -- стоимости ребра из стартовой вершины. Тогда давайте хранить в dist не точное расстояние d(start, v_i), а его оценку сверху. Теперь мы можем обновить значения в dist для всех вершин окружения.

for u, cu in g[start]:
    dist[u] = cu

Мы выставили значения каким-то вершинам. Что теперь? Посмотрим на нестартовую вершину v (v != start), имеющую наименьшую метку в dist. Учитывая, что отрицательных ребер в графе нет, то никакая из других вершин уже не сможет уменьшить оценку сверху для расстояния d(start, v), значит, расстояние d(start, v) совпадает со значением dist[v]. Однако ребра из вершины v вполне могут понизить оценки сверху для смежных с v вершин u_i до d(start, v) + c(v, u_i). Такое уменьшение оценки сверху на расстояние называется релаксацией.

for u, cu in g[v]:
    dist[u] = min(dist[u], dist[v] + cu)

После того, как мы рассмотрели вершину v, мы снова можем выбрать вершину v' с наименьшей оценкой сверху dist[v'], чтобы проводить дальнейшие релаксации. Однако нам нельзя выбирать уже просмотренные вершины start и v. Чтобы упростить эту проверку, заведем массив processed, изначально заполненный False, а выставлять значение True для вершины будем только в момент ее выбора как вершины с наименьшей оценкой сверху. Теперь процесс выбора следующей вершины заключается в поиске вершины с наименьшим значением в массиве dist и со значением False в массиве processed.

min_v = None
for v, dv in enumerate(dist):
    if not processed[v] and (min_v is None or dv < dist[min_v]):
        min_v = v

Когда наш алгоритм завершится? Можно заметить, что на каждой итерации мы ищем вершину из непросмотренных и выставляем ей метку, то есть ни одна вершина не может быть выбрана дважды. Это значит, что мы сделаем ровно n итераций, так как после уже не будет вершин, которые бы можно было обработать. После всех n итераций мы можем быть уверены, что массив dist хранит действительные расстояния до вершин.

Однако если граф несвязный, то в какой-то момент выбранной вершиной станет вершина с меткой INF в dist. Это значит, что дальше продолжать алгоритм нет смысла, ведь все последующие вершины, которые мы обработаем, не получат метки меньше INF, а значит, все они недостижимы из стартовой вершины.

if min_v is None or dist[min_v] == INF:
    break

Соберем все вместе:

def distances(start):
    dist[start] = 0
    
    for i in range(n):
        
        min_v = None
        for v, dv in enumerate(dist):
            if not processed[v] and (min_v is None or dv < dist[min_v]):
                min_v = v
        
        if min_v is None or dist[min_v] == INF:
            break
        
        processed[min_v] = True
        for u, cu in g[min_v]:
            dist[u] = min(dist[u], dist[min_v] + cu)

Можно заметить, что цикл for в данном случае излишний, так как в итерации n + 1 мы все равно не смогли бы найти следующую вершину и вышли бы из цикла, а переменную i мы нигде не используем. Избавимся от него:

def distances(start):
    dist[start] = 0
    
    while True: # <-- здесь!
        
        min_v = None
        for v, dv in enumerate(dist):
            if not processed[v] and (min_v is None or dv < dist[min_v]):
                min_v = v
        
        if min_v is None or dist[min_v] == INF:
            break
        
        processed[min_v] = True
        for u, cu in g[min_v]:
            dist[u] = min(dist[u], dist[min_v] + cu)

В программной реализации можно использовать и None вместо INF, но придется обрабатывать эти значения отдельно.

Время работы алгоритма составит O(n^2 + m), так как в худшем случае мы сделаем n итераций по O(n) каждая (просмотр массива dist) и просмотрим все ребра графа.

Более быстрый алгоритм

Какая часть в предыдущем алгоритме может быть оптимизирована? Мы не можем уменьшить число итераций, так как нам может понадобиться обойти все вершины. Мы не можем сократить число просмотренных ребер. Зато мы можем сделать поиск минимума быстрее.

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

def distances(start):
    q = PriorityQueue()
    q.enqueue((0, start))
    
    while not q.empty():
        (dv, v) = q.dequeue()
        
        if processed[v]:
            continue
        
        processed[v] = True
        dist[v] = dv
        
        for (u, cu) in g[v]:
            q.enqueue((dv + cu, u))

Так как мы в худшем случае просмотрим все ребра и для каждого ребра выполним операцию добавления элемента в очередь, то время работы алгоритма будет O(m log m).

Этот алгоритм называется алгоритмом Дейкстры.

Если мы хотим в дальнейшем восстановить кратчайшие пути, то для этого заведем массив pred, который будем обновлять при извлечении вершины из очереди. Для этого нам придется запоминать вершину-предка, из которой мы пришли, добавляя ее в очередь вместе с расстоянием как третье значение.

def distances(start):
    q = PriorityQueue()
    
    q.enqueue((0, start, None)) # <-- здесь!
    
    while not q.empty():
        (dv, v, pv) = q.dequeue()
        
        if processed[v]:
            continue
        
        processed[v] = True
        dist[v] = dv
        
        pred[v] = pv # <-- и здесь!
        
        for (u, cu) in g[v]:
        
            q.enqueue((dv + cu, u, v)) # <-- и здесь!

Взгляд с другой стороны

Вспомним, как выглядел алгоритм нахождения расстояний в ребрах с помощью BFS:

def distances(start):
    q = Queue()

    dist[start] = 0
    visited[start] = True
    q.enqueue(start)

    while not q.empty():
        v = q.dequeue()
        
        for u in g[v]:
            if not visited[u]:
                visited[u] = True
                dist[u] = dist[v] + 1
                q.enqueue(u)

Здесь мы обновляем значение dist для каждой вершины только один раз перед добавлением ее в очередь. Если мы добавляем вершину только один раз, то и извлекаем ее лишь один раз. Это значит, что мы можем обновлять массив dist и при извлечении вершины из очереди, но для этого нам нужно передавать само значение расстояния через очередь тоже:

def distances(start):
    q = Queue()

    visited[start] = True
    q.enqueue((0, start)) # <-- здесь!

    while not q.empty():
        dv, v = q.dequeue() # <-- и
        dist[v] = dv        # <-- здесь!
        
        for u in g[v]:
            if not visited[u]:
                visited[u] = True
                
                q.enqueue((dv + 1, u)) # <-- и еще здесь!

Также мы каждый раз проверяем вершину на посещенность перед добавлением в очередь. Но что будет, если ее не проверять? Вершина может попасть в очередь несколько раз. Какое же тогда значение будет верным? Чем позже вершина попала в очередь, тем больше расстояние, поэтому верным значением будет самое первое из значений, с которыми вершина попала в очередь. Это значит, что при извлечении вершины из очереди мы проверим, извлекали ли мы ее до этого, и если да, то просто проигнорируем. Посмотрим, что получилось:

def distances(start):
    q = Queue()
    
    q.enqueue((0, start)) # <-- убрали код здесь!

    while not q.empty():
        dv, v = q.dequeue()
        
        if visited[v]:      # <-- новый
            continue        # <-- блок
        visited[v] = True   # <-- здесь!
        
        dist[v] = dv
        
        for u in g[v]:                
            q.enqueue((dv + 1, u)) # <-- еще убрали здесь!

Теперь код выглядит аккуратнее, мы присваиваем метки только в одном месте (раньше для вершины start необходимо было присваивать метки отдельно), а также уровней вложенности стало меньше.

Сейчас вершины мы скорее не посещаем, а обрабатываем, поэтому заменим названием массива visited на processed:

def distances(start):
    q = Queue()
    q.enqueue((0, start))

    while not q.empty():
        dv, v = q.dequeue()
        if processed[v]:    # <-- здесь!
            continue
        
        processed[v] = True # <-- и здесь!
        dist[v] = dv
        
        for u in g[v]:                
            q.enqueue((dv + 1, u))

Посмотрим, что же находится внутри очереди. Очевидно, что ключи (расстояния до вершин) внутри очереди образуют неубывающую последовательность. Так как мы извлекаем элементы всегда из начала очереди, то, получается, мы всегда извлекаем минимум. Таким свойством обычная очередь обладает лишь засчет того, что мы добавляем элементы лишь на 1 больше извлеченного последним, в результате чего первый и последний элемент очереди отличаются не более, чем на 1. Однако мы знаем и более общую структуру, позволяющую нам выполнять те же операции: добавление элемента и извлечение минимума -- это приоритетная очередь. Давайте заменим обычную очередь на приоритетную:

def distances(start):
    q = PriorityQueue() # <-- здесь!
    q.enqueue((0, start))

    while not q.empty():
        dv, v = q.dequeue()
        if processed[v]:
            continue
        
        processed[v] = True
        dist[v] = dv
        
        for u in g[v]:                
            q.enqueue((dv + 1, u))

Заметим, что при добавлении вершины в очередь мы указываем новое расстояние как dv + 1, то есть расстояние до родителя, увеличенное на единицу. Эту единицу вполне можно рассматривать как вес ребра, а весь граф -- как взвешенный граф, все ребра которого имеют вес 1. А если в графе нет такого ограничения и вес ребер может быть любым неотрицательным числом? В этом случае новое расстояние будет вычисляться как dv + c(v, u), и нарушается предыдущее правило, что каждое последующее добавление увеличивает расстояние до вершины. Однако это уже не мешает алгоритму, и его работа ничуть не изменится: мы все еще будем помещать вершину в очередь несколько раз, а первое извлеченное расстояние для каждой вершины будет минимальным. Заменим невзвешенный граф на взвешенный:

def distances(start):
    q = PriorityQueue()
    q.enqueue((0, start))

    while not q.empty():
        dv, v = q.dequeue()
        if processed[v]:
            continue
        
        processed[v] = True
        dist[v] = dv
        
        for u, cu in g[v]:          # <-- вот
            q.enqueue((dv + cu, u)) # <-- здесь!

Сравните этот код с кодом из предыдущего пункта.

Отрицательные веса

До сих пор мы исходили из весьма оптимистичного допущения, что веса не могут быть отрицательными. Но почему это мешает алгоритму Дейкстры? Этот алгоритм основан на предположении, что если извлекли какое-то значение из очереди, то все дальнейшие добавленные в очередь значения будут не меньше этого, что очевидно неправда в случае с отрицательными весами.

Но что делать, если такая задача все же возникла? Казалось бы, очевидным решением будет просто прибавить ко всем весам какую-то достаточно большую константу, чтобы они стали неотрицательными, и запустить на новом графе Дейкстру. Да, в новом графе абсолютная разница между весами ребер сохранится, но мы ведь ищем кратчайшие пути, а абсолютная разница между их длинами изменится, ведь теперь она зависит и от их длин в ребрах. По этой причине такое решение может привести к неправильному ответу. В качестве примера можно взять граф из пяти вершин, который описывается двумя цепями a -- (-2) --> b -- (3) --> c и a -- (-2) --> d -- (1) --> e -- (1) --> c. В этом графе, очевидно, более коротким путем от a до c будет второй путь с расстоянием 0 (против 1 у первого пути), но если мы сделаем все веса неотрицательными, прибавив к ним константу 2, то кратчайшим путем станет первый путь с расстоянием 5 (против 6 у второго пути).

Для решения задачи воспользуемся еще одним допущением. Будем считать, что все кратчайшие пути содержат не более n - 1 ребер, так как ни один из них не содержит внутри цикла, то есть ни одна вершина не может встретиться в кратчайшем пути дважды. Если бы обратное было возможно, то в графе бы имелся цикл отрицательной стоимости (сумма весов ребер которого отрицательна), достижимый из вершины start, что позволило бы делать пути, проходящие через входящие в этот цикл вершины, сколь угодно маленькими за счет нужного числа проходов по этому циклу.

Вспомним наивную реализацию алгоритма Дейкстры и будем опять хранить в массиве dist оценки сверху на расстояние от стартовой вершины, изначально заполнив его INF. Тогда, так как кратчайшие пути имеют длину в ребрах не более n - 1, достаточно пройти n - 1 раз по всем ребрам, проводя релаксации, чтобы массив dist хранил действительные расстояния (с каждым проходом по ребрам мы "проталкиваем" правильный ответ все дальше и дальше от стартовой вершины). Если мы представим граф в виде списка ребер (v, c(v, u), u), то получим следующий алгоритм:

def distances(start):
    dist[start] = 0
    
    for i in range(n - 1):
        for v, c, u in edges:
            dist[u] = min(dist[u], dist[v] + c)

Если мы хотим в дальнейшем воостанавливать пути, то добавим массив pred, изначально заполненный None:

def distances(start):
    dist[start] = 0
    
    for i in range(n - 1):
        for v, c, u in edges:
            if dist[u] > dist[v] + c:  # <-- раскрыли
                dist[u] = dist[v] + c  # <-- код в
                pred[u] = v            # <-- блок `if`

Но при такой реализации мы так и не узнаем, верным ли было наше допущение об отсутствии отрицательных циклов или нет. Если оно было неверным, то и ответ мы получили неверный, поэтому желательно об этом сообщить. Как это понять? Давайте сделаем еще один проход по ребрам, и если где-то произошла релаксация, то какой-то кратчайший путь содержит как минимум n ребер, что говорит о наличии в нем цикла.

def distances(start):
    dist[start] = 0
    
    for i in range(n - 1):
        for v, c, u in edges:
            if dist[u] > dist[v] + c:
                dist[u] = dist[v] + c
                pred[u] = v
    
    for v, c, u in edges:          # <-- проверка на
        if dist[u] > dist[v] + c:  # <-- наличие
            panic!                 # <-- отрицательного цикла

Также, если на какой-либо итерации не было произведено ни одной итерации, то алгоритм можно остановить, ведь все последующие итерации никаких изменений тоже не внесут.

В итоге мы делаем n итераций по m ребрам и получаем время работы O(nm). Этот алгоритм называется алгоритмом Беллмана-Форда.

Если граф задан не списком ребер, а списком смежности, как ранее, то придется лишь поменять способ прохода по ребрам:

def distances(start):
    dist[start] = 0
    
    for i in range(n - 1):
        
        for v, gv in enumerate(g): # <-- вот
            for u, c in gv:        # <-- здесь!
                
                if dist[u] > dist[v] + c:
                    dist[u] = dist[v] + c
                    pred[u] = v
    
    for v, gv in enumerate(g): # <-- и вот
        for u, c in gv:        # <-- здесь!
            
            if dist[u] > dist[v] + c:
                panic!