k本の経路を求める簡単なアルゴリズムを LightGraphs.jl を用いて実装する

5 min read読了の目安(約4700字

背景

前回 LightGraphs.jl とその拡張である SimpleWeightedGraphs.jl について書いていました.

https://zenn.dev/takilog/articles/3012e77e242666649369

ここで多少は LightGraphs.jl の操作に慣れてきたので,例えば Dijkstra 法を使うには

g = # something graph
ds = dijkstra_shortest_paths(G, s)
distance = ds.dists[t]
path = enumerate_paths(ds)[t]

のような形で,スッと利用できると頭に染み込みました (これの使い方についてはドキュメントの Paths and Traversal を見て頂ければ良いです.あまりちゃんと書いてないですが…).

この記事ではこれらの基本的な機能を用いて,2点間のいくつかの経路を求めるアルゴリズムを実装してみました.最も古典的には k-shortest path problems (KSP) と呼ばれていますが,どのような k 本を選ぶかによって様々なバリエーションがあります.例えば東京駅から新宿駅まで移動する際にいろんな経路から好きな地下鉄を選ぶように,二点間を結ぶ経路をいくつか用意したいことは (よく) あります.

Yen's k-shortest path problems

k本の最短経路を求めるよく知られているアルゴリズムです.なんとLightGraphs.jlに実装されています!ドキュメントをチラ見してみましょう.

yen_k_shortest_paths(g, source, target, distmx=weights(g), K=1; maxdist=Inf);

Perform Yen's algorithm on a graph, computing k-shortest distances between source and target other vertices. Return a YenState that contains distances and paths.

お,使えるじゃん!と思ったところですが,残念ながらバグっていて使えないみたいです (使える場合もあるかもしれないですが,確かにこのバグは確認しました(2021/4現在)).

https://github.com/JuliaGraphs/LightGraphs.jl/issues/1505

仕方ないので (実はこのバグは後で気づきました?),自分で変わりの手法を実装します.

手法

ここでは2つの手法を実装してみます.例題のグラフはこのようなやつです.始点s=1と終点t=8を用い,通常の最短経路を描いています.辺の移動コストはそのまま平面上のユークリッド距離が設定されています.

経路ペナルティ

さて,Yen's algorithmや他のKSPのアルゴリズムを真面目に実装しても良いのですが,古典的なKSPのアルゴリズムは似た経路がたくさん含まれるような集合が返されることが知られています.ちょっと考えるとピンと来るのですが,最短経路をちょっと変えただけの経路はそこそこ短いコストになるので,それが解にはいってしまいます.これを避けるために真面目な研究がされていますが,ここでは雑な実装をやってみます.

ある経路 P_{s,t} を選択したとき,この P_{s,t} に既に含まれた辺 e\in P_{s,t} については,他の解には含まれなくても良さそうな気がします (本当でしょうか?).これを既存の Dijkstra 法をなるべく生かしたまま実装するには,e\in P_{s,t} の重み w_e を,パス P_{s, t} を計算する度にワザと大きく変化させれば良い気がします (結果としてこの辺を通ることが少なくなる気がします).式で書くと,適当な値 p をつかって長さを長くするように

w(e) \to w(e) \times (1 + p)

変化させれば良さそうです.これを実装しましょう.前回まで使っていた SimpleWeightedGraphs は辺自体に重み情報が入っていてスッとしていますが,残念ながらデータ構造としては immutable なので変更できません.そのため,一度グラフ G (これは通常のLightGraphs.jlでのグラフです) と復元した重み情報 \{w_e\} を行列で持ったものを用意し,行列では重みを操作し, Dijkstraで計算するときに変化させた重みを参照するようにします. 求める本数は K=5,ペナルティ倍率は p=0.5 とします.

s = 1
t = 8
p = 0.5
K = 5
dists = []
for k in 1:K
    # dmat情報でs-t最短経路を求める
    ds = dijkstra_shortest_paths(G, s, dmat)
    path_st = enumerate_paths(ds)[t]
    push!(paths, path_st)
    
    # ペナルティ
    for i in 1:length(path_st) - 1
        u, v = path_st[i], path_st[i + 1]
        nw = dmat[u, v] * (1 + p)
        dmat[u, v] = dmat[v, u] = nw
    end
end

ものすごく簡単な (というか何もしていない) 実装ですが,試した結果を示します.そこそこ動作している気がしますね.

1 (st) 2 3 4 5

グラフ重みランダム化

先に述べた手法では,計算した経路 P_{s,t} の上の情報を編集することで,既に計算した経路を避けるような次の経路を求めていました.また別のアプローチとして,一度グラフ上の最短経路を計算したら,グラフの辺情報をランダムに摂動させるという更に大胆なアプローチを考えます.全体がフルフル変化すれば,Dijkstraを適用して得られる経路も変化するだろうという真っ直ぐなお気持ちです.どのような変化をさせるかはいろいろありえますが,一番簡単に現在の重み w(e) を平均とした正規分布程度のぶれを発生 (ただしDijkstraを利用する場合には重みが負にならないように気をつける) させると良さそうです.

具体的には重みが負 (とついでに0) にならないようなパラメータ \tau > 0 と適当な値 \delta を用いて

w(e) \to \max\left(\tau, w(e) + \mathcal{N}(0, w(e)^2 \delta^2) \right)

のように変化させることを考えます.実装はほとんど同じで,経路ペナルティを計算していた箇所を

for e in edges(G)
    nw = max(τ, e.weight + e.weight ^ 2 * delta ^ 2 * randn())
    dmat[e.src, e.dst] = dmat[e.dst, e.src] = nw
end

のように重みを書き換えます.残りの処理は経路ペナルティ法と同じです.こちらについても結果を観察してみましょう.

1 2 3 4 5

肉眼ではほとんど差がわからない気もしますが,LightGraphs.jl の機能を十分に活かした状態で,できるだけ簡単に複数の経路を得ることができました.

参考情報