それでは毛玉諸君、これにて失敬

日々の精進を備忘録的に綴ります。

ベルマンフォード法とダイクストラ法の概念を完全に理解する

概要

社内の勉強会でダイクストラ法(dijkstra)とベルマンフォード法(bellman-ford)の紹介をしました。その際にアルゴリズムの証明を色々調べていましたが、解説されている記事があまり多くないことに気づきました。 「なぜそうなるの?」という部分は個人的に大事な要素だと思うので、その辺も踏まえてこれらのアルゴリズムを紹介したいと思います。

筆者はどんな人?

プログラミング歴:1年ちょい(Pythonのみ) Atcoder茶(よわよわ) 2020年4月からデータ分析委託の企業に勤めているが、なぜかアルゴリズムの勉強ばっかりやっている。

このアルゴリズムはいつ使うの?

重み付きグラフの単一始点最短経路問題を解く際に使います。

それぞれの用語ですが、 重み付き:辺を通るのに必要なコスト 単一始点:スタート地点が一か所 最短経路:スタート地点から目的地まで最小コストで向かう経路 です。

下のグラフを例にすると、 image.png スタート地点$A$から全ての辺への最短経路を求めることが目的になります。 ちなみに$A→E$の最短経路は $(A→D→B→F→E)$ です。これをアルゴリズムを使用して求めていきます。

アルゴリズムの前に、グラフの大事な性質

1.最短路の部分道は最短路

最短路の部分道はどこをとっても最短路になります。

<証明> image.png 仮に上のグラフで$A→E$の最短路を $(A→B→C→D→E cost:5)$ とします。

このうち$(A→B→C cost:3)$ はこの最短路の部分道ですが、この区間にはB'という頂点を利用すれば コストを1少なく移動できます。 するとこの区間の最短路は $(A→B'→C cost:2)$ となるため、$A→E$を $(A→B'→C→D→E cost:4)$ で移動できてしまい、仮定と矛盾します。 最短路の部分道は必ず最短路になります。

2.経路緩和性

始点に関連する頂点から順番に最短路を求めていくと、最終的にすべての最短経路が求められます。 これは結局、始点から徐々に求めた最短路がその先に到達する頂点の部分道となるためです。

ベルマンフォード法とは

全頂点間の最短路を求めます。ダイクストラ法と大きく異なる点は、負の重みが許容されることです。

最短路を求めるプロセスは、

  1. 全ての辺に注目し、その辺を通った先の頂点のコストを小さくできれば値を更新
  2.  1.のプロセスを(頂点の数-1)回行う です。

下の図だと、

image.png

1.のプロセスで、$A$を出発し、$C$にコスト$5$で到達できるので、$C:5$とメモします。 次に、$A→B$は$-1$で到達できるので$B:-1$、$B→C$にはコスト$-1+2=1$で到達できます。 この時、$B$を経由するとより小さなコストでCに到達できることが分かったので、$C:1$と更新します。 この作業を全ての辺に対して行います。これが1回目。 あとは(頂点の数-1)回繰り返せば、$A$から出発して全ての頂点への最短路を求めることが出来ます。

なぜ(頂点の数-1)回行うの?

始点から目的地までに通る頂点の最大の個数は、始点を除く頂点、すなわち(頂点の数-1)個です。 前述したプロセスでは、1回目のループで始点から移動できる頂点)の最短路が必ず求められます。 2回目のループでは(1回目で決まったルートから移動できる頂点)の最短路が必ず求められます。 ・ ・ ・ 結局このループは、始点をスタートして順番に最短路を求める作業になるので、経路緩和性により最短路を求めることができます。

負の重みによって困る場合

下の図を見てください。 image.png 丸で囲んだ部分はサイクルになっていますが、ここを1周するとコストが-2下がります。 つまりこのサイクルを使えば無限にコストを下げることが可能です。 これを負閉路といいます。

ベルマンフォード法では、負閉路の検出も可能です。 先ほど、(頂点の数-1)回のループで最短路を求めることを確認しましたが、再度全ての頂点のコストを計算し、コストを更新できる頂点がある場合は負閉路が存在するということになります。

いいからソースコード見せろよ

隣接リスト形式(競プロとかでよく利用する、隣接行列というのもある)の入力に対応し、ベルマンフォード法により最短路を検出します。 コードの下に負閉路のある入力例を示します(上の図のグラフ)。負閉路を通る頂点のコストを$-inf$に置き換えています。

クリックしたら見れるよ

def shortest_path(s,n,es):
    #s→iの最短距離
    # s:始点, n:頂点数, es[i]: [辺の始点,辺の終点,辺のコスト]
    
    #重みの初期化。ここに最短距離をメモしていく
    d = [float("inf")] * n
    
    #スタート地点は0
    d[s] = 0

    #最短路を計算
    for _ in range(n-1):
        # es: 辺iについて [from,to,cost]
        for p,q,r in es:
            #矢印の根が探索済、かつコストができるルートがあれば更新
            if d[p] != float("inf") and d[q] > d[p] + r:
                d[q] = d[p] + r
    
    #負閉路をチェック
    for _ in range(n-1):
        # e: 辺iについて [from,to,cost]
        for p,q,r in es:
            #更新される点は負閉路の影響を受けるので-infを入れる
            if d[q] > d[p] + r:
                d[q] = -float("inf")
    return d
################
#入力
n,w = map(int,input().split()) #n:頂点数 w:辺の数

es = [] #es[i]: [辺の始点,辺の終点,辺のコスト]
for _ in range(w):
    x,y,z = map(int,input().split())
    es.append([x,y,z])
    # es.append([y,x,z])

print(shortest_path(0,n,es))

'''
入力
7 12
0 1 2 
0 2 -5
0 3 5
1 2 3
1 4 1
1 5 3
2 4 6
3 2 -2
3 5 4
4 6 -5
5 4 -2
6 2 -3
'''
>>[0, 2, -inf, 5, -inf, 5, -inf]

計算量は$O(V2+E)$(頂点の数:$V$、辺の数:$E$)です。

ダイクストラ法とは

非負の重みを持つグラフで使用できます。

最短路を求めるプロセスはヨビノリさんの動画 )が死ぬほど分かりやすいです。 (ヨビノリさんの守備範囲広すぎる…)

一応こちらでもご説明します。 image.png

頂点のコストを(初期値:$S = 0$, その他$ = ∞$)とします。

  1. 一番コストが小さい頂点を始点とします。最初は$S$が始点です。
  2. 始点のコストを確定します(図では緑になる)。この頂点のコストはこれ以降更新されません。
  3. 始点から移動可能な頂点に出発点のコスト+辺重みを頂点にメモします。中段列の頂点に$7、3、1$とメモされました。コストはまだ更新されるかもしれないので未確定です。

  4. すべての頂点が確定するまで1~3を繰り返します。メモされた頂点の中で1のコストを持つ頂点が一番小さいので確定し、その頂点から移動可能な頂点をメモします。もし、より小さなコストで移動できるのならコストを更新します。

なぜ最短路が求められる?

未確定の頂点のうち最小のコストを持つ頂点は、それ以上小さなコストに更新することができないためです。 非負の重みしか存在しないので、別ルートを通るにしても0以上のコストが発生します。そのため未確定のうち、最小のコストは最短路と確定することが出来るのです。

で、ソースコードは?

クリックしたら見れるよ

import heapq

def dijkstra(s, n, es):
    #始点sから各頂点への最短距離
    #n:頂点数, cost[u][v] : 辺uvのコスト(存在しないときはinf)
    d = [float("inf")] * n
    #探索リスト
    S = [s]
    #優先度付きキューを使うと計算量が少なくなるよ
    heapq.heapify(S)

    d[s] = 0
    
    while S:
        #最小コストの頂点を取り出し
        v = heapq.heappop(S)
        # print("-----------v", v)
        for u, w in es[v]:
            # print(u, w)
            #より小さなコストで到達可能であれば更新し、リストに追加
            if d[u] > d[v]+w:
                d[u] = d[v]+w
                heapq.heappush(S, u)
        # print("d", d)
    return d

n,w = map(int,input().split()) #n:頂点数 w:辺の数

cost = [[float("inf") for i in range(n)] for i in range(n)] 
es = [[] for _ in range(n)]

for i in range(w):
    x, y, z = map(int, input().split())
    es[x].append([y, z])
    es[y].append([x, z]) #有向グラフの場合はこの行を消す
# print(es)
print(dijkstra(0, n, es))

オリジナルの計算量はベルマンフォード法と同じ$O(V2)$ですが、優先度付きキューを使用すれば$O((E+V)logV)$に削減できます。

参考文献とかサイト

アルゴリズムイントロダクション 第3版 第2巻: 高度な設計と解析手法・高度なデータ構造・グラフアルゴリズム (世界標準MIT教科書) (T.コルメン,R.リベスト,C.シュタイン,C.ライザーソン,浅野哲夫,岩野和生,梅尾博司,山下雅史,和田幸一) かなり重い本ですが、グラフ理論の証明がとても丁寧でした。

グラフ理論⑤(ダイクストラアルゴリズム) (予備校のノリで学ぶ「大学の数学・物理」) https://www.youtube.com/watch?v=X1AsMlJdiok こちらの動画を見るだけで、ダイクストラ法は完全に理解できます!

・蟻本 python 単一最短経路法2(ダイクストラ法) 競技プログラミング (じゅっぴーさんのブログ、コードはこちらを参考にさせて頂きました) https://juppy.hatenablog.com/entry/2018/11/01/%E8%9F%BB%E6%9C%AC_python%E5%8D%98%E4%B8%80%E6%9C%80%E7%9F%AD%E7%B5%8C%E8%B7%AF%E6%B3%952%EF%BC%88%E3%83%80%E3%82%A4%E3%82%AF%E3%82%B9%E3%83%88%E3%83%A9%E6%B3%95%EF%BC%89%E7%AB%B6%E6%8A%80%E3%83%97

おまけ

記事中のグラフはnetworkxで作図しました。 気になる人向けにコードを置いておきます。

クリックしたら見れるよ

#ダイクストラ法のグラフ
import matplotlib.pyplot as plt
import networkx as nx

#有向グラフのセット
G = nx.DiGraph()

#辺を追加していく
# G.add_edges_from([("A", "B", {"weight":2)])
# nx.add_path(G, ["A", "C"])
# G.edges["A", "B"]["weight"] = 2
G.add_weighted_edges_from([("A", "B", 5), ("A", "C", 7), ("A", "D", 1), ("B", "C", 3), ("B", "F", 3), ("D", "B", 2)])
G.add_weighted_edges_from([("C", "E", 3), ("D", "F", 8), ("F", "E", 2), ("B", "E", 6)])

#図中で頂点の座標を指定しないとランダムに配置される
pos1 = {}
pos1["A"] = (0, 0)
pos1["B"] = (1, 0)
pos1["C"] = (1, 2)
pos1["D"] = (1, -2)
pos1["E"] = (2, 1)
pos1["F"] = (2, -1)

#そのまま描写するとweightの文字が邪魔なので消す
edge_labels = {(i, j): w['weight'] for i, j, w in G.edges(data=True)}

#重みを書き込む
nx.draw_networkx_edge_labels(G, pos1, edge_labels=edge_labels, font_size=18)
#作図
nx.draw_networkx(G, pos1, node_size=1000, label_pos=100, node_color="r", font_size=18)
#図の保存
plt.savefig("dijkstra.png")
plt.show()

#networkxにもダイクストラ法が実装されている。計算量はO(V^2+E)なので遅い
pred, dist = nx.dijkstra_predecessor_and_distance(G, "A")
print(pred, dist)

>>{'A': [], 'B': ['D'], 'C': ['B'], 'D': ['A'], 'F': ['B'], 'E': ['F']} 
>>{'A': 0, 'D': 1, 'B': 3, 'C': 6, 'F': 6, 'E': 8}