SciPyで最小全域木を求める: minimum_spanning_tree
SciPyの関数scipy.sparse.csgraph.minimum_spanning_tree()を使うとグラフの最小全域木を取得できる。
ここでは以下の内容について説明する。
- 最小全域木と
scipy.sparse.csgraph.minimum_spanning_tree() minimum_spanning_tree()への入力- リスト、
numpy.ndarray scipy.sparse行列（疎行列）
- リスト、
- 出力された最小全域木グラフの処理
- 隣接行列への変換
- 重み（コスト）の総和
- 辺のリスト化
- 入力の隣接行列の扱い
- 処理速度比較（引数
overwriteの効果）
scipy.sparse行列（疎行列）についての詳細は以下の記事を参照。
以下の説明内容およびサンプルコードはSciPy1.3.1のもの。バージョンが違うと挙動が異なる場合があるので注意。
なお、scipy.sparse.csgraph.minimum_spanning_tree()はSciPy0.11.0で追加された。AtCoderでも使える。
最小全域木とscipy.sparse.csgraph.minimum_spanning_tree()
ある無向グラフの部分グラフで、すべての頂点を含む木を全域木という。木は閉路を持たない連結グラフ（任意の2頂点間に経路が存在するグラフ）。
さらに、重み付きグラフにおいて辺の重み（コスト）の総和が最小となる全域木を最小全域木という
例えば下図の左のグラフに対して、右の木が最小全域木となる。
input graph minimum spanning tree
(0) (0)
/ \ /
3 8 3
/ \ /
(3)---5---(1) (3)---5---(1)
\ / /
6 2 2
\ / /
(2) (2)
上図の左のグラフに対して、下図の左は全域木であるが最小全域木でない（重みが最小でない）。右は木でない（連結されていない頂点が存在する）。
input graph minimum spanning tree
(0) (0)
/ \ /
3 8 3
/ \ /
(3)---5---(1) (3)---5---(1)
\ / /
6 2 2
\ / /
(2) (2)
SciPyの関数scipy.sparse.csgraph.minimum_spanning_tree()を使うとグラフから最小全域木を求めることができる。最初の図は以下のページから引用した。
最小全域木を求めるアルゴリズムにはプリム法とクラスカル法があるが、minimum_spanning_tree()ではクラスカル法を使っている。クラスカル法の計算量は、頂点の数をV、辺の数をEとしたとき、O(E log V)またはO(E log E)となる。
アルゴリズムや計算量についての詳細は以下のWikipediaのページを参照。
minimum_spanning_tree()への入力
minimum_spanning_tree()の第一引数にグラフの隣接行列を指定する。2次元リスト（リストのリスト）やNumPy配列ndarray、scipy.sparse行列（疎行列）を指定可能。
なお、いずれの場合も内部でscipy.sparse.csr_matrixに変換されるので、最初からcsr_matrixを生成できるのであればそちらのほうが効率的。
以下のようにインポートしておく。
import numpy as np
from scipy.sparse.csgraph import minimum_spanning_tree
from scipy.sparse import csr_matrix, coo_matrix, lil_matrix
リスト、numpy.ndarray
上で図示したグラフの隣接行列は以下のようになる。
l = [[0, 8, 0, 3],
[0, 0, 2, 5],
[0, 0, 0, 6],
[0, 0, 0, 0]]
これをminimum_spanning_tree()の第一引数に指定すると最小全域木が得られる。返り値はscipy.sparse.csr_matrix。
mst = minimum_spanning_tree(l)
print(mst)
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
print(type(mst))
# <class 'scipy.sparse.csr.csr_matrix'>
結果のcsr_matrixの処理については後述。
numpy.ndarrayでもOK。
a = np.array(l)
print(a)
# [[0 8 0 3]
# [0 0 2 5]
# [0 0 0 6]
# [0 0 0 0]]
print(type(a))
# <class 'numpy.ndarray'>
print(minimum_spanning_tree(a))
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
scipy.sparse行列（疎行列）
scipy.sparseにはcsr_matrixやcoo_matrixなど様々な格納形式の疎行列のクラスがある。これらを引数に指定することもできる。
csr = csr_matrix(l)
print(csr)
# (0, 1) 8
# (0, 3) 3
# (1, 2) 2
# (1, 3) 5
# (2, 3) 6
print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>
print(minimum_spanning_tree(csr))
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
どの形式でもOK。上述のように、内部でcsr_matrixに変換される。
print(minimum_spanning_tree(coo_matrix(l)))
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
print(minimum_spanning_tree(lil_matrix(l)))
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
例えばcsr_matrixオブジェクトは、辺の重み、および、その辺を構成する頂点のリストから生成できる。
n = 4
d = [8, 3, 2, 5, 6]
i = [0, 0, 1, 1, 2]
j = [1, 3, 2, 3, 3]
csr_ = csr_matrix((d, (i, j)), shape=(n, n))
print(csr_)
# (0, 1) 8
# (0, 3) 3
# (1, 2) 2
# (1, 3) 5
# (2, 3) 6
print(csr_.toarray())
# [[0 8 0 3]
# [0 0 2 5]
# [0 0 0 6]
# [0 0 0 0]]
print(minimum_spanning_tree(csr))
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
元のグラフのデータがこのような辺の2つの頂点とその重みで与えられた場合は、リストやnumpy.ndarrayで2次元配列を生成するよりもcsr_matrixを生成するほうが楽。
scipy.sparse行列の生成についての詳細は以下の記事を参照。
出力された最小全域木グラフの処理
minimum_spanning_tree()は結果をcsr_matrixで返す。
mst = minimum_spanning_tree(l)
print(mst)
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
print(type(mst))
# <class 'scipy.sparse.csr.csr_matrix'>
隣接行列への変換
csr_matrixのtoarray()でnumpy.ndarrayに変換できる。
print(mst.toarray())
# [[0. 0. 0. 3.]
# [0. 0. 2. 5.]
# [0. 0. 0. 0.]
# [0. 0. 0. 0.]]
値は浮動小数点数floatなので、整数intに変換したい場合はastype()を使う。csr_matrixにもnumpy.ndarrayにもastype()がある。
print(mst.toarray().astype(int))
# [[0 0 0 3]
# [0 0 2 5]
# [0 0 0 0]
# [0 0 0 0]]
print(mst.astype(int).toarray())
# [[0 0 0 3]
# [0 0 2 5]
# [0 0 0 0]
# [0 0 0 0]]
リストに変換したい場合はさらにnumpy.ndarrayのtolist()を使う。
print(mst.toarray().astype(int).tolist())
# [[0, 0, 0, 3], [0, 0, 2, 5], [0, 0, 0, 0], [0, 0, 0, 0]]
print(type(mst.toarray().astype(int).tolist()))
# <class 'list'>
重み（コスト）の総和
重み（コスト）の総和はcsr_matrixのsum()メソッドを使う。
print(mst.sum())
# 10.0
print(int(mst.sum()))
# 10
辺のリスト化
csr_matrixのnonzero()で最小全域木の辺を構成する頂点のペアでリスト化できる。
r, c = mst.nonzero()
print(r, c)
# [0 1 1] [3 2 3]
print(list(zip(*mst.nonzero())))
# [(0, 3), (1, 2), (1, 3)]
重みはdata属性に格納されている。
print(mst.data)
# [3. 2. 5.]
print(list(zip(*mst.nonzero(), mst.data.astype(int))))
# [(0, 3, 3), (1, 2, 2), (1, 3, 5)]
組み込み関数zip()と*を使った転置については以下の記事を参照。
入力の隣接行列の扱い
minimum_spanning_tree()は入出力のグラフを無向グラフとして扱う。隣接行列の[i, j]と[j, i]がいずれも非0の場合、小さい方が辺の重みとして使われる。
l = [[0, 8, 0, 3],
[8, 0, 2, 5],
[0, 2, 0, 6],
[3, 5, 6, 0]]
print(minimum_spanning_tree(l).toarray().astype(int))
# [[0 0 0 3]
# [0 0 2 5]
# [0 0 0 0]
# [0 0 0 0]]
l = [[0, 8, 0, 3],
[100, 0, 2, 5],
[100, 100, 0, 6],
[100, 100, 100, 0]]
print(minimum_spanning_tree(l).toarray().astype(int))
# [[0 0 0 3]
# [0 0 2 5]
# [0 0 0 0]
# [0 0 0 0]]
l = [[0, 8, 0, 3],
[0, 0, 100, 5],
[0, 2, 0, 6],
[0, 0, 0, 0]]
print(minimum_spanning_tree(l).toarray().astype(int))
# [[0 0 0 3]
# [0 0 0 5]
# [0 2 0 0]
# [0 0 0 0]]
値が同じ場合は上の方（行番号が小さい方）が使われる模様。
処理速度比較（引数overwriteの効果）
minimum_spanning_tree()の第二引数overwriteはTrueにすると効率的というような説明がなされている。デフォルトはFalse。
overwrite : bool, optional
if true, then parts of the input graph will be overwritten for efficiency.
scipy.sparse.csgraph.minimum_spanning_tree — SciPy v1.3.0 Reference Guide
これの効果を確認する。
結論を先に書いておくと、第一引数への入力がcsr_matrixであるときのみ第二引数overwriteは効果がある。すべての場合で高速になるわけではない。また、データ入力時のムダなコピーを減らすだけなので、処理時間が大幅に削減されるわけでもない。
以下の例はJupyter Notebookのマジックコマンド%%timeitを利用しており、Pythonスクリプトとして実行しても計測されない。
引数overwriteに関連するソースコードは以下。
第一引数に入力された値がscipy.sparse行列でない場合はoverwriteは処理に使われていない。したがって、リストやnumpy.ndarrayの場合はoverwriteをTrueにしても変わらない。
l = [[0, 8, 0, 3],
[0, 0, 2, 5],
[0, 0, 0, 6],
[0, 0, 0, 0]]
%%timeit
minimum_spanning_tree(l)
# 378 µs ± 9.78 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
minimum_spanning_tree(l, True)
# 383 µs ± 16.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
第一引数に入力された値がscipy.sparse行列であったとき、csr_matrixへの変換を行うcsr_matrix()の引数copyにnot overwriteが渡されている。overwrite=Trueのときcopy=False、overwrite=False（デフォルト）のときcopy=True。
csr_matrix()の引数copy=Falseが効くのは入力がcsr_matrixのときのみ。したがって、第一引数に入力された値がcsr_matrixでなければoverwriteをTrueにしても変わらない。
csr_matrixだと速くなる。
csr = csr_matrix(l)
%%timeit
minimum_spanning_tree(csr)
# 158 µs ± 5.95 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%%timeit
minimum_spanning_tree(csr, True)
# 108 µs ± 2.19 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
coo = coo_matrix(l)
%%timeit
minimum_spanning_tree(coo)
# 185 µs ± 11.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
minimum_spanning_tree(coo, True)
# 184 µs ± 9.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
overwrite=Trueとすると元の行列の値が上書きされてしまいそうな印象だがそのようなことはない。最小全域木自体は別の行列として生成される。
mst = minimum_spanning_tree(csr, True)
print(csr)
# (0, 3) 8
# (1, 2) 3
# (1, 3) 2
# (1, 3) 5
# (2, 3) 6
print(mst)
# (0, 3) 3.0
# (1, 2) 2.0
# (1, 3) 5.0
元がリストやnumpy.ndarrayの場合、csr_matrixに変換してからoverwrite=Trueとしてminimum_spanning_tree()を実行すると速くなる。
%%timeit
minimum_spanning_tree(csr_matrix(l))
# 364 µs ± 4.11 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
minimum_spanning_tree(csr_matrix(l), True)
# 319 µs ± 11.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
元のグラフのデータがどういった形で提供されるかによるが、リストやnumpy.ndarrayを経由せずにcsr_matrixを生成する手段があればそちらのほうがよいだろう。
なお、グラフのサイズが大きくなると最小全域木を求める処理自体に時間がかかるため、overwriteの効果は相対的に小さくなっていくものと思われる。