note.nkmk.me

SciPyでグラフの最短経路を算出(ダイクストラ、ベルマンフォードなど)

Date: 2019-08-25 / tags: Python, SciPy

scipy.sparse.csgraphの関数shortest_path()を使うとグラフの最短経路問題を解くことができる。単一始点最短経路問題にも全点対最短経路問題にも対応。

ここでは以下の内容について説明する。

  • 最短経路問題
  • 各アルゴリズムに対応: shortest_path()
    • 基本的な使い方
    • 単一始点最短経路問題: 引数indices
    • アルゴリズムの指定: 引数method
    • 有向グラフか無向グラフか: 引数directed
    • 重み付きか重みなしか: 引数unweighted
    • 経路復元: 引数return_predecessors
  • ワーシャルフロイド法: floyd_warshall()
  • ダイクストラ法: dijkstra()
  • ベルマンフォード法: bellman_ford()
  • ジョンソン法: johnson()
  • 処理速度比較
    • nupmy.ndarrayscipy.sparse
    • shortest_path()と個別の関数
    • 引数indicesの効果
    • ベルマンフォード法とジョンソン法

shortest_path()はデフォルトで入力されたグラフに応じて適切なアルゴリズムを自動で選択してくれる。使用するアルゴリズムが決まっている場合は、floyd_warshall()dijkstra()などの個別の関数を実行したほうが無駄なバリデーションなどを省略できるので高速。

引数などの使い方はshortest_path()も個別のアルゴリズムの関数も同じ。

以下の説明内容およびサンプルコードはSciPy1.3.1のもの。バージョンが違うと異なる場合があるので注意。

以下のようにインポートする。

import numpy as np
from scipy.sparse.csgraph import shortest_path, floyd_warshall, dijkstra, bellman_ford, johnson
from scipy.sparse import csr_matrix
スポンサーリンク

最短経路問題

最短経路問題(最短路問題)の説明は以下の通り。

グラフ理論における最短経路問題(さいたんけいろもんだい、英: shortest path problem)とは、重み付きグラフの与えられた2つのノード間を結ぶ経路の中で、重みが最小の経路を求める最適化問題である。
最短経路問題 - Wikipedia

特定の頂点(ノード)とその他の頂点との最短経路問題を「単一始点最短経路問題」、任意の2頂点の組み合わせの最短経路問題を「全点対最短経路問題」という。

各アルゴリズムに対応: shortest_path()

はじめに、各アルゴリズムに対応したshortest_path()について説明する。

基本的な使い方

以下の隣接行列のグラフを例とする。

l = [[0, 1, 2, 0],
     [0, 0, 0, 1],
     [3, 0, 0, 3],
     [0, 0, 0, 0]]
 ┌---(0)<--┐
 1       ↓2│3↑
 ↓         ↓
(1)       (2)
 │         │
 1         3
 └-->(3)<--┘

shortest_path()の第一引数に指定する。リストだとエラー。NumPy配列numpy.ndarrayはOK。返り値はnumpy.ndarrayとなる。

# print(shortest_path(l))
# AttributeError: 'list' object has no attribute 'shape'

a = np.array(l)
print(type(a))
# <class 'numpy.ndarray'>

print(shortest_path(a))
# [[ 0.  1.  2.  2.]
#  [inf  0. inf  1.]
#  [ 3.  4.  0.  3.]
#  [inf inf inf  0.]]

print(type(shortest_path(a)))
# <class 'numpy.ndarray'>

デフォルトでは全点対最短経路問題に対する解が返される。各要素の値が最短経路のコストの総和(最短距離)となっている。

例の場合、0から1への最短経路のコストの総和は[0, 1]の要素の値12から1への最短経路のコストの総和は[2, 1]の要素の値4。経路が存在しない頂点の組み合わせの値はinfとなる。

最短経路のパス(道のり)を取得するには引数return_predecessorsTrueにして返り値を処理する。後述。

shortest_path()の第一引数にはnumpy.ndarrayのほかscipy.sparse行列(疎行列)も指定できる。

csr = csr_matrix(l)
print(csr)
#   (0, 1)  1
#   (0, 2)  2
#   (1, 3)  1
#   (2, 0)  3
#   (2, 3)  3

print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>

print(shortest_path(csr))
# [[ 0.  1.  2.  2.]
#  [inf  0. inf  1.]
#  [ 3.  4.  0.  3.]
#  [inf inf inf  0.]]

scipy.sparse行列の生成についての詳細は以下の記事を参照。

最後に紹介するが、csr_matrixを指定したほうがnumpy.ndarrayよりも高速。

単一始点最短経路問題: 引数indices

上述のように、デフォルトでは全点対最短経路問題に対する解となり、任意の2頂点の組み合わせの最短距離が求められる。

単一始点最短経路問題を解きたい、つまり、始点(または終点)が決まっている場合は、引数indicesを指定する。

整数intまたはそのリストで指定する。

print(shortest_path(csr, indices=0))
# [0. 1. 2. 2.]

print(shortest_path(csr, indices=[0, 3]))
# [[ 0.  1.  2.  2.]
#  [inf inf inf  0.]]

アルゴリズムの指定: 引数method

最短経路問題を解くアルゴリズムはいくつかある。shortest_path()ではデフォルトで第一引数に入力されたグラフから適切なアルゴリズムを選択して実行する。

引数methodで任意のアルゴリズムを指定することも可能。

  • 'auto'(デフォルト)
    • 入力グラフに応じて自動的に選択
  • 'FW'
    • ワーシャルフロイド法(Floyd-Warshall algorithm)
    • O[N^3]
  • 'D'
    • ダイクストラ法(Dijkstra’s algorithm)
      • フィボナッチヒープ使用
    • O[N(N*k + N*log(N))]
    • 負の重みがあると使えない
  • 'BF'
    • ベルマンフォード法(Bellman-Ford algorithm)
    • O[N(N^2 k)]
  • 'J'
    • ジョンソン法(Johnson’s algorithm)
    • ベルマンフォード法とダイクストラ法の組み合わせ

計算量は公式ドキュメントに記載の値で、kは頂点に接続される辺の平均数。全点対最短経路問題時のものであることに注意。

'auto'の場合の判定基準はソースコードを参照。

以下のように選択される。

  • 引数indicesが指定されている、または、第一引数の非ゼロ要素数が全要素数の25%より少ない
    • 負の重みがある
      • → ジョンソン法
    • 負の重みがない
      • → ダイクストラ法
  • それ以外 (= 引数indicesが指定されている、かつ、第一引数の非ゼロ要素数が全要素数の25%以上)
    • → ワーシャルフロイド法

なお、使いたいアルゴリズムが決まっている場合、shortest_path()method='D'などとするよりもdijkstra()などの関数を直接実行したほうが無駄な処理が省略できて高速。最後に実験結果を紹介する。

有向グラフか無向グラフか: 引数directed

引数directed=Falseとすると無向グラフとして処理される。デフォルトはTrue(有向グラフ)。

print(shortest_path(csr, directed=False))
# [[0. 1. 2. 2.]
#  [1. 0. 3. 1.]
#  [2. 3. 0. 3.]
#  [2. 1. 3. 0.]]

例の場合、以下のようなグラフと等価とみなされる。

l_ud = [[0, 1, 2, 0],
        [1, 0, 0, 1],
        [2, 0, 0, 3],
        [0, 1, 3, 0]]

print(shortest_path(csr_matrix(l_ud)))
# [[0. 1. 2. 2.]
#  [1. 0. 3. 1.]
#  [2. 3. 0. 3.]
#  [2. 1. 3. 0.]]

directed=Falseではグラフの[i, j][j, i]がどちらも非ゼロ要素の場合、小さい方の値がその頂点間の距離となる。

l2 = [[0, 1, 2, 0],
      [100, 0, 0, 1],
      [100, 0, 0, 3],
      [0, 100, 100, 0]]

print(shortest_path(csr_matrix(l2), directed=False))
# [[0. 1. 2. 2.]
#  [1. 0. 3. 1.]
#  [2. 3. 0. 3.]
#  [2. 1. 3. 0.]]

ダイクストラ法とジョンソン法ではdirected=False[i, j][j, i]が等価でない場合に結果がおかしくなる可能性があるそうなので注意。

As currently implemented, Dijkstra’s algorithm and Johnson’s algorithm do not work for graphs with direction-dependent distances when directed == False. i.e., if csgraph[i,j] and csgraph[j,i] are non-equal edges, method=’D’ may yield an incorrect result.
scipy.sparse.csgraph.shortest_path — SciPy v1.3.0 Reference Guide

重み付きか重みなしか: 引数unweighted

引数unweighted=Trueとすると重みなしのグラフとして処理される。デフォルトはFalse(重み付き)。

print(shortest_path(csr, unweighted=True))
# [[ 0.  1.  1.  2.]
#  [inf  0. inf  1.]
#  [ 1.  2.  0.  1.]
#  [inf inf inf  0.]]

すべての重みが1となることと等価。

l_uw = [[0, 1, 1, 0],
        [0, 0, 0, 1],
        [1, 0, 0, 1],
        [0, 0, 0, 0]]

print(shortest_path(csr_matrix(l_uw)))
# [[ 0.  1.  1.  2.]
#  [inf  0. inf  1.]
#  [ 1.  2.  0.  1.]
#  [inf inf inf  0.]]

経路復元: 引数return_predecessors

最短経路のコストの総和(最短距離)だけでなく最短経路自体を取得したい場合、引数return_predecessorsTrueとする。

以下のように、これまでの結果と合わせて、最短経路において一つ前の頂点を示すpredecessorsnumpy.ndarrayで返される。

d, p = shortest_path(csr, return_predecessors=True)

print(d)
# [[ 0.  1.  2.  2.]
#  [inf  0. inf  1.]
#  [ 3.  4.  0.  3.]
#  [inf inf inf  0.]]

print(p)
# [[-9999     0     0     1]
#  [-9999 -9999 -9999     1]
#  [    2     0 -9999     2]
#  [-9999 -9999 -9999 -9999]]

predecessorsを辿っていくことで最短経路を取得できる。

例えば以下のような関数を定義する。もっと効率的な方法があるかもしれない。

def get_path(start, goal, pred):
    return get_path_row(start, goal, pred[start])

def get_path_row(start, goal, pred_row):
    path = []
    i = goal
    while i != start and i >= 0:
        path.append(i)
        i = pred_row[i]
    if i < 0:
        return []
    path.append(i)
    return path[::-1]

結果は以下の通り。経路がない場合は空のリストが返される。

print(get_path(0, 3, p))
# [0, 1, 3]

print(get_path(2, 1, p))
# [2, 0, 1]

print(get_path(3, 3, p))
# [3]

print(get_path(3, 1, p))
# []

引数indicesを指定した場合はpredecessorsもその頂点のものだけになる。

d2, p2 = shortest_path(csr, indices=2, return_predecessors=True)

print(d2)
# [3. 4. 0. 3.]

print(p2)
# [    2     0 -9999     2]

print(get_path_row(2, 1, p2))
# [2, 0, 1]

print(get_path_row(2, 0, p2))
# [2, 0]

ワーシャルフロイド法: floyd_warshall()

ここから個別のアルゴリズムの関数について説明する。

ワーシャルフロイド法(Warshall–Floyd Algorithm)は全点対最短経路問題を解くアルゴリズム。

全点対最短経路問題のためのものなので引数indicesは指定できない。

# print(floyd_warshall(csr, indices=0))
# TypeError: floyd_warshall() got an unexpected keyword argument 'indices'

# print(shortest_path(csr, indices=0, method='FW'))
# ValueError: Cannot specify indices with method == 'FW'.

ダイクストラ法: dijkstra()

ダイクストラ法(Dijkstra's algorithm)は辺の重みが非負数のグラフの単一始点最短経路問題を解くアルゴリズム。

scipy.sparse.csgraphdijkstra()ではフィボナッチヒープを利用して実装されている。

辺の重みに負の値が含まれている場合は警告が出る。この例ではコメントアウトしているが、エラーではなく警告であるため実行できてしまう。後述のように負の閉路がある場合は特に注意。

l_n = [[0, 1, 2, 0],
       [0, 0, 0, 1],
       [-3, 0, 0, 3],
       [0, 0, 0, 0]]

# print(dijkstra(csr_matrix(l_n)))
# UserWarning: Graph has negative weights: dijkstra will give inaccurate results
#              if the graph contains negative cycles. Consider johnson or bellman_ford.

shortest_path()にはない引数limitが使用可能(バージョン0.14.0以降)。limitの値よりコストの総和が大きくなると探索を打ち切る。

print(dijkstra(csr))
# [[ 0.  1.  2.  2.]
#  [inf  0. inf  1.]
#  [ 3.  4.  0.  3.]
#  [inf inf inf  0.]]

print(dijkstra(csr, limit=2))
# [[ 0.  1.  2.  2.]
#  [inf  0. inf  1.]
#  [inf inf  0. inf]
#  [inf inf inf  0.]]

ベルマンフォード法: bellman_ford()

ベルマンフォード法(Bellman–Ford algorithm)はグラフの単一始点最短経路問題を解くアルゴリズム。辺の重みは負数でもよい。

グラフに負の閉路(negative cycle)が存在すると最短経路を求めることができない。

ベルマンフォード法では負の閉路を検出してエラーNegativeCycleErrorを出す。

l_nc = [[0, 1, 2, 0],
        [0, 0, -10, 1],
        [3, 0, 0, 3],
        [0, 0, 0, 0]]

# print(bellman_ford(csr_matrix(l_nc)))
# NegativeCycleError: Negative cycle detected on node 0

ワーシャルフロイド法、ジョンソン法でもエラーとなるが、ダイクストラ法では負の閉路を検出できないため警告が出るのみで正しくない結果が返されてしまう。

# print(floyd_warshall(csr_matrix(l_nc)))
# NegativeCycleError: Negative cycle in nodes [0 1 2]

# print(johnson(csr_matrix(l_nc)))
# NegativeCycleError: Negative cycle detected on node 0

print(dijkstra(csr_matrix(l_nc)))
# [[  0.   1.  -9.  -6.]
#  [ -7.   0. -10.  -7.]
#  [  3.   4.   0.   3.]
#  [ inf  inf  inf   0.]]
# 
# /usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: Graph has negative weights: dijkstra will give inaccurate results if the graph contains negative cycles. Consider johnson or bellman_ford.
#   """Entry point for launching an IPython kernel.

例外処理を行う場合は以下の記事を参照。

ジョンソン法

ジョンソン法(Dijkstra's algorithm)はグラフの単一始点最短経路問題を解くアルゴリズム。辺の重みは負数でもよい

ダイクストラ法とベルマンフォード法を組み合わせたもの(らしい)。

辺の重みが非負数であればダイクストラ法のほうが速いが、ジョンソン法は負の閉路を検出できるというメリットがある。

処理速度比較

いくつかの場合の処理速度(処理時間)を比較する。

以下のようなグラフを作成して例とする。100個の頂点に対して約200個の辺をランダムに生成している。辺の数がぴったり200個でないのはランダム生成の結果、同じ行・列に値が存在したから。

import numpy as np
from scipy.sparse.csgraph import shortest_path, dijkstra, floyd_warshall, bellman_ford, johnson
from scipy.sparse import csr_matrix

n = 100
c = n * 2
np.random.seed(1)
d = np.random.randint(0, n, c)
i = np.random.randint(0, n, c)
j = np.random.randint(0, n, c)

csr = csr_matrix((d, (i, j)), shape=(n, n))
a = csr.toarray()

print(a.shape)
# (100, 100)

print((a > 0).sum())
# 198

結果としては、第一引数をscipy.sparse.csr_matrixとして、shortest_path()ではなく個別の関数を呼ぶのが速い。

なお、当然ながら頂点および辺の個数やその構成によって結果は異なるのであくまでも参考程度とされたい。

以下の例はJupyter Notebookのマジックコマンド%%timeitを利用しており、Pythonスクリプトとして実行しても計測されないので注意。

nupmy.ndarrayとscipy.sparse

内部での変換処理が省略されるため、scipy.sparse.csr_matrixを第一引数としたほうがnumpy.ndarrayより速い。

%%timeit
shortest_path(a)
# 2.3 ms ± 73.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit
shortest_path(csr)
# 1.52 ms ± 21.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

ソースコードをざっと見た限りfloyd_warshall()では内部処理をnumpy.ndarrayで行っているようだが、後述のように、floyd_warshall()でもcsr_matrixのほうが速い。何か勘違いをしているかもしれない。

shortest_path()と個別の関数

上の例では、引数indicesを指定せず、第一引数の非ゼロ要素数が全要素数の25%より少ないのでダイクストラ法が選択される。

shortest_path()methodで指定した場合は自動選択の場合とほぼ同じ。

%%timeit
shortest_path(a, method='D')
# 2.25 ms ± 65.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit
shortest_path(csr, method='D')
# 1.51 ms ± 25.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

個別の関数dijkstra()のほうが若干速い。

%%timeit
dijkstra(a)
# 1.9 ms ± 52.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit
dijkstra(csr)
# 1.4 ms ± 13.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

なお、この例の場合はダイクストラ法よりワーシャルフロイド法のほうが速い。shortest_path()のアルゴリズム自動選択が常に最良の結果とは限らない。

%%timeit
shortest_path(a, method='FW')
# 1.11 ms ± 33.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%%timeit
shortest_path(csr, method='FW')
# 555 µs ± 71.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%%timeit
floyd_warshall(a)
# 737 µs ± 13.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%%timeit
floyd_warshall(csr)
# 433 µs ± 22.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

引数indicesの効果

始点(または終点)が決まっている場合は引数indicesを指定したほうが当然ながら速い。

%%timeit
dijkstra(a, indices=0)
# 590 µs ± 28.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%%timeit
dijkstra(csr, indices=0)
# 111 µs ± 6.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

ただし、この例では頂点100個から引数indicesで1個だけを選択しているが、処理時間が1 / 100になるわけではない。

特に頂点が少ないグラフ(サイズの小さいグラフ)の場合、メインの探索処理のほかの前処理などにかかる時間が相対的に長くなるため、indicesで結果を絞っても大幅に処理時間が削減されるということにはならない。

ベルマンフォード法とジョンソン法

辺の重みに負の値を含む場合。

ジョンソン法のほうがベルマンフォード法よりも速いが、この例の場合はワーシャルフロイド法のほうがさらに速い。

a_n = a.copy()
a_n[0, 1] = -10
csr_n = csr_matrix(a_n)

%%timeit
johnson(csr_n)
# 1.5 ms ± 17.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%%timeit
bellman_ford(csr_n)
# 6.12 ms ± 141 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit
floyd_warshall(csr_n)
# 402 µs ± 16.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

上の例と同じく、始点(または終点)が決まっている場合は引数indicesを指定したほうが速い。

%%timeit
johnson(csr_n, indices=0)
# 201 µs ± 6.84 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
スポンサーリンク
シェア
このエントリーをはてなブックマークに追加

関連カテゴリー

関連記事