全組最短路

合抱の木も毫末より生じ、九層の台も累土より起こり、千里の行も足下より始まる。

――老子, 『老子道徳経』 64 章 (紀元前六世紀)

And I would walk five hundred miles,

And I would walk five hundred more,

Just to be the man who walks a thousand miles

To fall down at your door.

[五百マイル歩き、

また五百マイル歩く。

ただ千マイル歩いて、

君の家の前で倒れる男になるために。]

――The Proclaimers, "I'm Gonna Be (500 Miles)", Sunshine on Leith (2001)

もう少し...もう少し...

――Red Leader [Drewe Henley], Star Wars (1977)

導入

前章では単一の頂点 \(s\) から他の全ての頂点への最短路を \(s\) を根とした最短路木を作ることで見つけるアルゴリズムについて議論しました。最短路木はグラフの各頂点 \(v\) に対して二つの情報を割り当てます:

この章で扱うのは、もっと一般的な全組最短路問題 (all pairs shortest path problem, APSP) です。これは任意のソース頂点から任意の目標頂点への最短路を求める問題であり、この問題を解くアルゴリズムは \(G\) の全ての頂点の組 \(u, v\) に対して次の情報を計算します:

この直感的な定義は二つの例外的なケースをいくつか無視していますが、どれも前章で見たものです:

全組最短路問題の出力は二つの \(V \times V\) の配列であり、一つには\(V^{2}\) 個の最短路の距離2が、もう一つには \(V^{2}\) 個の前者が格納されます。ですがこの章では、距離の配列を計算することに集中します。最短路を実際に構築するために必要になる前者配列は、アルゴリズムをほんの少し改変するだけで計算できます (ヒント、ヒント)。

たくさんの単一ソース

全組最短路問題の明らかな解法は、単一ソース最短路問題を各頂点に対して一回ずつ、合計で \(V\) 回解くというものです。具体的に言うと、単一ソース最短路問題を解くたびに一次元の部分配列 \(\mathit{dist}[s, \cdot]\) が埋まります:

procedure \(\texttt{ObviousAPSP}\)(\(V, E, w\))

for 頂点 \(s\) do

\(\mathit{dist}[s, \cdot] \leftarrow \)\(\texttt{SSSP}\)(\(V, E, w, s\))

このアルゴリズムの実行時間が用いる単一ソース最短路アルゴリズムに依存しているのも明らかです。単一ソースのときと同じように、グラフの構造と辺の重みに依存する四種類の選択肢があります:

重みの変更

負の辺があると実行が非常に遅くなります: どうにかして削除できないでしょうか? たくさんの人が思い浮かべる単純なアイデアは、全ての辺の重みを同じだけ増加させることで全ての辺の重みを正にして、Bellman-Ford の代わりに Dijkstra を使うというものです。残念ながら、この単純なアイデアは上手く行きません ――こうすると、多くの辺を持つ路の重みが少ない辺しか持たない路よりも多く増加してしまうからです。次の図に示すように、全ての辺を同じ比率で増加させると、辺をたくさん持つ路は辺を少ししか持たない路よりも速いペースで大きくなります。そのため二つの頂点の間の最短路が変化してしまう可能性があります。

全ての辺の重みを 2 増やすと \(s\) から \(t\) への最短路が変化する
図 9.1
全ての辺の重みを 2 増やすと \(s\) から \(t\) への最短路が変化する

しかしもっと手の込んだ方法を使うと、最短路を保存したまま重みを変更できます。この重みの変更方法の発見者はこの方法を応用した最短路アルゴリズムを 1973 年に示した Donald Johnson (ドナルド・ジョンソン) とされることが多いですが、実は Johnson はこの方法を 1972 年に発表された Jack Edmonds (ジャック・エドモンズ) と Richard Karp (リチャード・カープ) による論文から取ったとしています。また Nobuaki Tomizawa (冨澤 信明) も 1971 年に同じ手法を発表しているほか、Delbert Fulkerson (デルバート・ファルカーソン) による 1961 年の論文でも同じ手法が少し違った形式で示されています。

各頂点 \(v\) に対して、価値 (price) \(\pi (v)\) が割り振られているとします。頂点の価値は正でも負でもゼロでもあり得ます。このとき、新しい重み関数 \(w^{\prime}\) を次のように定義します: \[ w^{\prime}(u \rightarrow v) = \pi(u) + w(u \rightarrow v) - \pi (v) \] 直感的な説明をすると、\(u\) を離れるときには出国税 \(\pi (u)\) を払い、\(v\) に入るときには入国ギフトとして \(\pi(v)\) もらえるということです。

新しい重み関数 \(w^{\prime}\) に対する最短路が最初の重み関数 \(w\) に対する最短路と同じであることを示すのは難しくありません。実際に証明するとこうなります: 任意の路 \(u \rightsquigarrow v \) の全ての中間頂点 \(x\) について、入国ギフトとして \(\pi(x)\) をもらえますが、すぐに \(x\) から出ていくので、このギフトは出国税 \(\pi(x)\) で打ち消されます。したがって \[ w^{\prime}(u \rightsquigarrow v) = \pi(u) + w(u \rightsquigarrow v) - \pi (v) \] が成り立ちます。\(u\) から \(v\) への全ての路の長さが同じ値だけずれるので、\(u\) から \(v\) への最短路は変化しません (違う頂点の組を結ぶ路の長さのずれは異なるので、そのような路同士の長さの上下関係は変わる可能性があります)。

Johnson のアルゴリズム

Johnson (ジョンソン) の全組最短路アルゴリズムは全ての辺の重みが非負となるように各頂点に価値 \(\pi(v)\) を割り当て、その後 Dijkstra のアルゴリズムを使って新しい重みに対する最短路を計算します。

まずグラフの全ての頂点に到達できる頂点 \(s\) が存在すると仮定します。Johnson のアルゴリズムは最初 \(s\) から他の全ての頂点に対する最短路を (負辺があっても正しく動く) Bellman-Ford のアルゴリズムを使って計算し、それから価値関数を \(\pi(v) = \mathit{dist}(s, v)\) と定義します。こうした上で、新しい辺の重みを \[ w^{\prime}(u \rightarrow v) = \mathit{dist}(s, u) + w(u \rightarrow v) - \mathit{dist}(s, v) \] と定義します。Bellman-Ford が停止しているので、この新しい辺の重みは非負です! 辺 \(u \rightarrow v\) が緊張しているとは \(\mathit{dist}(s, u) + w(u \rightarrow v) < \mathit{dist}(s, v)\) が成り立つことであり、単一ソース最短路アルゴリズムは緊張している辺を全て取り除くことを思い出してください (Bellman-Ford が負閉路を検出した場合には最短路が well-defined でないので、Johnson のアルゴリズムはエラーを出して終了します)。

全ての頂点に到達するという条件を満たす頂点 \(s\) が存在しない場合、Bellman-Ford をどこから始めたとしても、いくつかの頂点の価値が無限大になってしまいます。この問題を避けるために、Johnson のアルゴリズムはどんな場合でも最初にまず新しい頂点 \(s\) をグラフに追加して、\(s\) から他の全ての頂点に重さゼロの辺を追加します。このとき他の頂点から \(s\) に戻る辺は追加しません。このように \(s\) を追加すれば \(s\) が中間点となる路は生まれないので、全ての組の間の最短路は変化しません。

Johnson のアルゴリズムの完全な疑似コードを次に示します。アルゴリズムの実行時間は Dijkstra のアルゴリズムの呼び出しによって支配されます。具体的には、\(\textsc{BellmanFord}\) を呼ぶのに \(O(VE)\) 時間かかり、\(\textsc{Dijkstra}\) を \(V\) 回呼ぶのに \(O(VE\log V)\) かかり、他のデータの管理に \(O(V + E)\) 時間かかるので、全体の実行時間は \(\pmb{O(VE \log V)} = O(V^{3} \log V)\)4 となります。結局、負辺があっても遅くなりません!

procedure \(\texttt{JohnsonAPSP}\)(\(V, E, w\))

⟨⟨人工的な頂点を追加する⟩⟩

新しい頂点 \(s\) を追加する

for \(s\) 以外の頂点 \(v\) do

新しい辺 \(s \rightarrow v\) を追加する

\(w(s \rightarrow v) \leftarrow 0\)

⟨⟨頂点の価値を計算する⟩⟩

\(\mathit{dist}[s, \cdot] \leftarrow \)\(\texttt{BellmanFord}\)(\(V, E, w, s\))

if \(\textsc{BellmanFord}\) が負閉路を検出した then

エラーを出す

⟨⟨辺の重みを変更する⟩⟩

for \((u,v) \in E\) do

\(w^{\prime}(u \rightarrow v) \leftarrow \mathit{dist}[s, u] + w(u \rightarrow v) - \mathit{dist}[s, v]\)

⟨⟨重みの変更されたグラフに対する最短路を計算する⟩⟩

for 頂点\(u\) do

\(\mathit{dist}^{\prime}[u, \cdot] \leftarrow \)\(\texttt{Dijkstra}\)(\(V, E, w^{\prime}, u\))

⟨⟨元のグラフに対する最短路を計算する⟩⟩

for 頂点\(u\) do

for 頂点\(v\) do

\(\mathit{dist}[u,v] \leftarrow \mathit{dist}^{\prime}[u, v] - \mathit{dist}[s, u] + \mathit{dist}[s, v]\)

図 9.2
Johnson の全組最短路アルゴリズム

動的計画法

単一ソース最短路アルゴリズムの代わりに動的計画法を使っても全組最短路問題を解くことができます。\(E = \Theta(V^{2})\) であるような密なグラフに対しては、動的計画法を使ったアプローチの方が Johnson のアルゴリズムよりも単純で (少しだけ) 速いアルゴリズムとなります。 この章ではこれから、入力グラフに負閉路が無いことを仮定します。

まずは動的計画法のアルゴリズムに付き物の再帰方程式を考えます。単一ソースと同様に、"明らかな" 再帰的定義 \[ \mathit{dist}(u, v) = \begin{cases} 0 & \text{if } u = v \\ \min\limits_{x \rightarrow v} (\mathit{dist}(u, x) + w(x \rightarrow v)) & \text{それ以外} \end{cases} \] は入力グラフが DAG であるときにしか正しくありません。有向閉路があると無限ループになってしまうためです。

そこで Bellman-Ford のアルゴリズムと同じように、追加のパラメータを追加してこの無限ループを停止させます。\(\mathit{dist}(u, v, l)\) で \(u\) から \(v\) への辺が \(l\) 個以下の歩道の中で最短のものの長さを表すことにします。任意の二つの頂点の間の最短路は最大でも \(V-1\) 個の辺しか持たないことから、真の最短路の長さは \(\mathit{dist}(u, v, V-1)\) です。単一ソース経路問題に対する Bellman の再帰方程式はこの関数にも適用できます: \[ \mathit{dist}(u, v, l) = \begin{cases} 0 & \text{if } l = 0 \text{ and } u = v \\ \infty & \text{if } l = 0 \text{ and } u \neq v \\ \min \left\lbrace \begin{array}{c} \mathit{dist}(u, v, l-1) \\ \min\limits_{x \rightarrow v} (\mathit{dist}(u, x, l-1) + w(x \rightarrow v)) \end{array} \right\rbrace & \text{ それ以外} \end{cases} \]

この再帰方程式を動的計画法のアルゴリズムに変形するのは簡単です。このアルゴリズムの実行時間は \(\pmb{O(V^{2}E)} = O(V^{4})\) となります。

procedure \(\texttt{ShimbelAPSP}\)(\(V, E, w\))

for 頂点 \(u\) do

for 頂点 \(v\) do

if \(u = v\) then

\(\mathit{dist}[u, v, 0] \leftarrow 0\)

else

\(\mathit{dist}[u, v, 0] \leftarrow \infty \)

for \(l \leftarrow 1\) to \(V-1\) do

for 頂点 \(u\) do

for 頂点 \(v \neq u\) do

\(\mathit{dist}[u, v, l] \leftarrow \mathit{dist}[u, v, l - 1]\)

for \(x \rightarrow v\) の形をした辺 do

if \(\mathit{dist}[u, v, l] > \mathit{dist}[u, x, l-1] + w(x \rightarrow v)\) then

\(\mathit{dist}[u, v, l] \leftarrow \mathit{dist}[u, x, l-1] + w(x \rightarrow v)\)

このアルゴリズムは 1954 年に Alfonso Shimbel (アルフォンソ・シンベル) によって初めて示されました5。Bellman による Bellman-Ford のアルゴリズムの定式化で見たように、内側の \(l\) と \(v\) に対するループは取り除くことができます。実際に取り除いたアルゴリズムを次に示します:

procedure \(\texttt{AllPairsBellmanFord}\)(\(V, E, w\))

for 頂点 \(u\) do

for 頂点 \(v\) do

if \(u = v\) then

\(\mathit{dist}[u, v, 0] \leftarrow 0\)

else

\(\mathit{dist}[u, v, 0] \leftarrow \infty \)

for \(l \leftarrow 1\) to \(V-1\) do

for 頂点 \(u\) do

for \(x \rightarrow v\) の形をした辺 do

if \(\mathit{dist}[u, v, l] > \mathit{dist}[u, x, l-1] + w(x \rightarrow v)\) then

\(\mathit{dist}[u, v, l] \leftarrow \mathit{dist}[u, x, l-1] + w(x \rightarrow v)\)

どのように導出したとしても、最終的に完成するアルゴリズムが \(V\) 個の異なる Bellman-Ford の実行をバラバラに並べたものと全く同じであるというのは驚くに値しません。具体的に言うと、メイン for ループの \(l\) 回目の反復が終わったとき、\(\mathit{dist}[u, v, l]\) は \(u\) から \(v\) への辺を \(l\) 個以下しか使わない路の中で最短のものの長さを表します。

分割統治

1971 年に Michael Fischer (マイケル・フィッシャー) と Albert Meyer (アルバート・マイヤー) によって提案された方法を使うと、さらに格段な高速化が可能です。Bellman の再帰方程式は最短路を少しだけ短い路と一つの辺に分けて考え、目標頂点の全ての前者を考えることで最短路を求めていました。こうする代わりに、最短路をその中間地点にある頂点で分けて二つの路にしてみましょう。こうすると、先ほど定義した \(\mathit{dist}(u, v, l)\) に対する異なる再帰方程式が手に入ります。中間の頂点が存在するのは \(l = 2\) までなので、ベースケースは \(l = 0\) ではなく \(l = 1\) です。また再帰方程式を簡単にするために、全ての頂点 \(v\) に対して \(w (v \rightarrow v) = 0\) と定義しておきます。 \[ \mathit{dist}(u, v, l) = \begin{cases} w(u \rightarrow v) & \text{if } l = 1 \\ \min\limits_{x} (\mathit{dist}(u, x, l/2) + \mathit{dist}(x, v, l/2)) & \text{それ以外} \end{cases} \]

この再帰方程式のままだと \(l\) が \(2\) のべきであるときにしか計算を行うことができず、それ以外の場合には (最大でも) 分数個の辺を持つ最短路を計算することになります。しかしこれは問題ではありません。なぜなら \(\mathit{dist}(u, v, l)\) は 任意の \(l \geq V - 1\) に対して \(u\) から \(v\) の最短路の長さと等しいからです。よって \(l = 2^{\lceil \lg V \rceil} < 2V\) とすれば上手く行きます。

ここでも、動的計画法のアルゴリズムはこの再帰方程式から簡単に手に入ります。実行時間が \(\pmb{O(V^{3} \log V)}\) であることは、アルゴリズムを実際に書かなくても式を見るだけで分かります ――\(u, v, x\) に対して \(V\) 個の値を考え、\(l\) に対しては \(\lg V\) 個の値を考えるからです。次に示す Fischer と Meyer のアルゴリズムの疑似コードでは、\(\mathit{dist}[u,v,i]\) という配列の要素が \(\mathit{dist}(u, v, 2^{i})\) の値を表します。

procedure \(\texttt{FischerMeyerAPSP}\)(\(V, E, w\))

for 頂点 \(u\) do

for 頂点 \(v\) do

\(\mathit{dist}[u, v, 0] \leftarrow w(u \rightarrow v)\)

⟨⟨\(l = 2^{i}\)⟩⟩

for \(i \leftarrow 1\) to \(\lceil \lg V \rceil\) do

for 頂点 \(u\) do

for 頂点 \(v\) do

\(\mathit{dist}[u, v, i] \leftarrow \infty\)

for 頂点 \(x\) do

if \(\mathit{dist}[u, v, i] > \mathit{dist}[u, x, i-1] + \mathit{dist}[x, v, i-1]\) then

\(\mathit{dist}[u,v,i] \leftarrow \mathit{dist}[u, x, i-1] + \mathit{dist}[x, v, i-1]\)

これまでのアルゴリズムが単一ソース最短路アルゴリズムを \(V\) 回実行するのと同じであったのに対して、このアルゴリズムはそうではありません。つまり、一番内側の処理は辺を緩和しているわけではありません。

Bellman-Ford を初めとするこれまでの動的計画法でそうしてきたように、このアルゴリズムでも配列の最後の次元を取り除くことができ、\(\mathit{dist}[u, v, i]\) の部分を \(\mathit{dist}[u, v]\) とすることで消費空間量を \(O(V^{3})\) から \(O(V^{2})\) に落とせます。この洗練されたアルゴリズムは 1957 年に Leyzorek et al. によって、彼らが Dijkstra のアルゴリズムを示したのと同じ論文で発表されました。

procedure \(\texttt{LeyzorekAPSP}\)(\(V, E, w\))

for 頂点 \(u\) do

for 頂点 \(v\) do

\(\mathit{dist}[u, v] \leftarrow w(u \rightarrow v)\)

\(\color{maroon}{\langle \langle l = 2^{i} \rangle \rangle}\)

for \(i \leftarrow 1\) to \(\lceil \lg V \rceil\) do

for 頂点 \(u\) do

for 頂点 \(v\) do

for 頂点 \(x\) do

if \(\mathit{dist}[u, v] > \mathit{dist}[u, x] + \mathit{dist}[x, v]\) then

\(\mathit{dist}[u, v] \leftarrow \mathit{dist}[u, x] + \mathit{dist}[x, v]\)

おかしな行列積

有向グラフの最短路の計算と正方行列のべき乗の計算の間には密接な関係があります。この関係は Alfonso によって最初に発見され、後に Bellman も独立に発見しました。\(n \times n\) の行列の二乗を求めるアルゴリズムを次に示します。これを \(\textsc{FischerMeyerAPSP}\) の内側のループと比べてみてください (似ていることが分かりやすいように二つ目のアルゴリズムの記法を少し変えています)。

procedure \(\texttt{MatrixSquare}\)(\(A\))

for \(i \leftarrow 1\) to \(n\) do

for \(j \leftarrow 1\) to \(n\) do

\(A^{\prime}[i, j] \leftarrow 0\)

for \(k \leftarrow 1\) to \(n\) do

\(A^{\prime}[i, j] \leftarrow A^{\prime}[i, j] + A[i, k] \cdot A[k, j]\)

procedure \(\texttt{FischerMeyerInnerLoop}\)(\(D\))

for 頂点 \(u\) do

for 頂点 \(v\) do

\(D^{\prime}[u, v] \leftarrow \infty\)

for 頂点 \(x\) do

\(D^{\prime}[u, v] \leftarrow \min \lbrace D^{\prime}[u, v],\ D[u, x] + D[x, v] \rbrace\)

二つのアルゴリズムの唯一の違いは、一つ目のアルゴリズムの乗算の部分で二つ目のアルゴリズムは加算を使っていて、一つ目のアルゴリズムの加算の部分で二つ目のアルゴリズムは最小化を使っている点です。このことから、最短路アルゴリズムの内側のループは「min-plus 行列積」や「距離行列積」あるいは「おかしな行列積 (funny matrix multiplication)」と呼ばれることがあります。

低速な方のアルゴリズム \(\textsc{ShimbelAPSP}\) は重み行列 \(w\) の \(V-1\) "min-plus 乗" を計算するための標準的な反復アルゴリズムです。最初のループで対角線上に \(0\) が、それ以外の場所に \(\infty\) が並んだ min-plus 単位行列が作られ、二つ目のループが次の "min-plus" 乗を計算します。そして \(\textsc{FischerMeyerAPSP}\) ではこの反復法を二乗の繰り返しに置き換えています。これはちょうど一章で紹介した乗算アルゴリズムと同じであり、エジプト人 \(\acute{\alpha} \rho \pi \varepsilon \delta o \nu \acute{\alpha} \pi \tau \alpha \iota\) の影響をここにも見ることができます!

整数乗算に対する Karatsuba の分割統治アルゴリズムと同じように、(通常の) 行列乗算に対するさらに高速な分割統治アルゴリズムが存在します。そのようなアルゴリズムの中で最初に示されたのは 1969 年の Volker Strassen (フォルカー・シュトラッセン) によるアルゴリズムであり、\(n \times n\) の行列乗算を七つの \(n/2 \times n/2\) の行列乗算で計算するものです。Strassen のアルゴリズムの実行時間は \(O(n^{\lg 7}) = O(n^{2.807355})\) です。Strassen のアルゴリズムは過去五十年間の間に何度も改善され、2018 年時点で最速の行列乗算アルゴリズムの実行時間は \(O(n^{2.372864})\) です。

二乗の繰り返しよりも高速なこれらのアルゴリズムは全て減算を使っているのですが、残念ながら "おかしな" 減算なるものは存在しません (\(\min\) の逆の演算とは何でしょうか?)。そのため少なくとも一般的なグラフに対しては、動的計画法を使ったアルゴリズムの内側のループを高速化する明らかな方法は存在しません。

しかし "明らかな方法がない" は "不可能である" を意味しません!実は、全組最短路問題の特殊ケースに対する格段に高速なアルゴリズムがいくつか存在します。そのようなアルゴリズムの中で最も素晴らしいものの一つは、1991 年に Zvi Galil と Oded Margalit によって発見され、1992 年に Raimund Seidel によって簡略化された単純な乱択アルゴリズムであり、このアルゴリズムを使うと重み無し無向グラフの全組最短路を平均 \(O(M(V) \log V)\) 時間で計算できます。ここで \(M(n) = O(n^{2.372864})\) は \(n \times n\) の整数行列の (通常の意味の) 乗算にかかる時間です6。Galil, Margalit, Seidel のアプローチはそれからさらに拡張され、今では小さい重みを持つ有向グラフの実際の最短路を、三乗を大きく下回る時間で決定的に計算できるようになっています。

重みが小さい整数の場合には大きな進展がみられているものの、辺に一般的な重みの付いたグラフに対する全組最短路問題が \(O(V^{2.999999})\) 時間で解けるのかどうかは、9 の数に関わらず誰も知りません。さらにそのようなアルゴリズムが不可能であることを示唆する証拠もいくつか存在します! そのためおそらくは、"明らかな方法がない" は "不可能である" を本当に意味するのでしょう。

(Kleene-Roy-) Floyd-Warshall (-Ingerman)

Johnson のアルゴリズムのフィボナッチヒープを使った実装と比べると、高速な方の動的計画法アルゴリズムでさえ最悪ケースの計算時間が \(O(\log V)\) 倍低速です。最短路問題の異なる定式化を使うとこの対数部分を改善できます。この定式化を使ったアルゴリズムは 1962 年に Robert Floyd (ピーター・フロイト) と Peter Ingerman (ピーター・インガーマン) によってぞれぞれ独立に発見され、両者とも同年の前半に Stephen Warshall (ステファン・ワーシャル) によって発表されたアルゴリズムを一般化したと記しています。また正確に言えば、この Warshall のアルゴリズムは 1959 年には Bernand Roy (バーナード・ロイ) によって、そして同じ再帰パターンは 1951 年に Stephen Kleene (スティーヴン・クリーネ)7 によって発見されています。

Warshall (と Roy と Kleene) は動的計画法の再帰方程式の三つ目のパラメータとして利用できる別の値を提案しました。この新しいパラメータは辺の数を制限するためではなく、通り抜けられる頂点を制限するために使われます。ここで「通り抜ける」とは「頂点に入ってかつ出る」を意味します。例えば \(w \rightarrow x \rightarrow y \rightarrow z\) という路は \(w\) から始まり、\(x, y\) を通り抜け、\(z\) で終わります。

頂点に \(1\) から \(V\) までの数字を適当に割り当てます。任意の頂点の組 \(u, v\) と任意の整数 \(r\) について、\(\pi(u, v, r)\) を次のように定義します:

\(u\) から \(v\) への路で通り抜ける頂点の数字が全て \(r\) 以下のものが存在するならば、そのような路の中で最短のものを \(\pi(u, v, r)\) とする。

この定義から \(\pi(u, v, V)\) が \(u\) から \(v\) への真の最短路となります。Kleene と Roy と Warshall は、この路についての単純な再帰的構造を発見しました:

制限された最短路 \(\pi(u, v, r)\) の再帰的構造
図 9.3
制限された最短路 \(\pi(u, v, r)\) の再帰的構造

\(\mathit{dist}(u, v, r)\) で \(\pi(u, v, r)\) の長さを表すことにすると、\(\pi(u, v, r)\) の再帰的な構造から直ちに次の再帰方程式が得られます: \[ \mathit{dist}(u, v, r) = \begin{cases} w(u \rightarrow v) & \text{if } r = 0 \\ \min \left\lbrace \begin{array}{c} \mathit{dist}(u, v, r-1) \\ \min (\mathit{dist}(u, r, r-1) + \mathit{dist}(r, v, r-1)) \end{array} \right\rbrace & \text{それ以外} \end{cases} \]

問題を解くために必要なのは、全ての頂点の組 \(u, v\) に対する \(\mathit{dist}(u, v, V)\) の計算です。この再帰方程式も動的計画法のアルゴリズムへと簡単に変換できます:

procedure \(\texttt{KleeneASPS}\)(\(V, E, w\))

for 頂点 \(u\) do

for 頂点 \(v\) do

\(\mathit{dist}[u,v,0] \leftarrow w(u \rightarrow v)\)

for \(r \leftarrow 1\) to \(V\) do

for 頂点 \(u\) do

for 頂点 \(v\) do

if \(\mathit{dist}[u,v,r-1] < \mathit{dist}[u,r,r-1] + \mathit{dist}[r,v,r-1]\) then

\(\mathit{dist}[u,v,r] \leftarrow \mathit{dist}[u,v,r-1]\)

else

\(\mathit{dist}[u,v,r] \leftarrow \mathit{dist}[u,r,r-1] + \mathit{dist}[r,v,r-1]\)

これまでに見てきた最短路問題に対する動的計画法のアルゴリズムと同じように、\(\textsc{KleeneASPS}\) もメモ配列の最後の次元を削除することで単純化できます。また頂点の番号付けは任意であることから、疑似コードにおいて番号を具体的に取得する必要はありません。これで、Floyd が改善した形の Warshall のアルゴリズムが手に入ります:

procedure \(\texttt{FloydWarshall}\)(\(V, E, w\))

for 頂点 \(u\) do

for 頂点 \(v\) do

\(\mathit{dist}[u,v] \leftarrow w(u \rightarrow v)\)

for 頂点 \(r\) do

for 頂点 \(u\) do

for 頂点 \(v\) do

if \(\mathit{dist}[u,v] < \mathit{dist}[u,r] + \mathit{dist}[r,v]\) then

\(\mathit{dist}[u,v] \leftarrow \mathit{dist}[u,r] + \mathit{dist}[r,v]\)

\(\textsc{FloydWarshall}\) と先述の少し遅い動的計画法アルゴリズム \(\textsc{LeyzorekAPSP}\) を比較すると興味深いことが分かります。\(\textsc{LeyzorekAPSP}\) では全ての頂点の三つ組を \(O(\log V)\) 回反復するのに対して \(\textsc{FloydWarshall}\) では一度の反復で済むのですが、二つのアルゴリズムの違いはループの順番だけなのです!

練習問題

    1. \(\textsc{LeyzorekAPSP}\) を改変して、最短路の長さの配列だけはなく前者へのポインタの配列も計算するようにしてください。実行時間は \(O(V^{3}\log V)\) のままであるようにしてください。

    2. \(\textsc{FloydWarshall}\) を改変して、最短路の長さの配列だけはなく前者へのポインタの配列も計算するようにしてください。実行時間は \(O(V^{3})\) のままであるようにしてください。

  1. この章で議論した全てのアルゴリズムは負閉路が含まれるグラフに対して正しい答えを計算できません。グラフに負閉路が含まれている場合、Johnson のアルゴリズムでは最初の Bellman-Ford が負閉路を検出して実行がそこで終了し、動的計画法のアルゴリズムでは間違った答えが計算されます。

    しかしこの章で議論した全てのアルゴリズムは、負閉路が含まれている場合でも最短路の長さを正しく計算するようにできます。つまり:

    • \(u\) が \(v\) に到達できないなら、\(\mathit{dist}[u, v] = \infty\) を返す。

    • \(u\) が \(v\) に到達できる負閉路に到達できるなら、\(\mathit{dist}[u, v] = -\infty\) を返す。

    • それ以外の場合、\(u\) から \(v\) への最短路が存在するので、その長さを返す。

    となるようにアルゴリズムを改変できるということです。

    1. \(\textsc{JohnsonAPSP}\) を改変して、負閉路が含まれている場合でも最短路の長さを正しく計算できるようにしてください。

    2. \(\textsc{LeyzorekAPSP}\) を改変して、負閉路が含まれている場合でも最短路の長さを正しく計算できるようにしてください。

    3. \(\textsc{FloydWarshall}\) を改変して、負閉路が含まれている場合でも最短路の長さを正しく計算できるようにしてください。

  2. この章で議論した全てのアルゴリズムは、入力グラフ \(G\) に負閉路が含まれているかどうかを報告するだけではなく、\(G\) に含まれる負閉路を検出したときにその負閉路を返すように改変できます。

    1. \(\textsc{JohnsonAPSP}\) を改変して、要素が最短路の長さの配列または負閉路を返すようにしてください。

    2. \(\textsc{LeyzorekAPSP}\) を改変して、要素が最短路の長さの配列または負閉路を返すようにしてください。

    3. \(\textsc{FloydWarshall}\) を改変して、要素が最短路の長さの配列または負閉路を返すようにしてください。

    全ての場合において、入力グラフに複数の負閉路が含まれる場合にはどれを選んでも構いません。

  3. \(G = (V, E)\) を辺に重みの付いた有向グラフとし、辺の重みは任意の実数で、負閉路は含まれないとします。

    1. 重みがゼロである閉路を見つけるか、そのような閉路が無ければそうだと正しく報告する効率の良いアルゴリズムを説明してください。

    2. 次の条件を満たす \(G\) の部分グラフ \(H\) を構築する効率の良いアルゴリズムを説明してください。

      • \(H\) は \(G\) の全ての頂点を含む。
      • \(H\) の全ての閉路は重みがゼロである。
      • \(G\) に含まれる重みがゼロの閉路は全て \(H\) に含まれる。

      例えば \(G\) がゼロ閉路を持たない場合、\(H\) は辺を持ちません。

  4. \(G = (V, E)\) を辺に重みの付いた有向グラフとし、辺の重みは任意の実数とします。\(G\) の頂点を \(k\) 個の互いに素な部分集合 \(V_{1}, V_{2}, \ldots, V_{k}\) に分割したとします。つまり、\(G\) の全ての頂点がちょうど一つの部分集合 \(V_{i}\) に属するということです。任意の添え字の組 \(i, j\) について、\(\delta(i, j)\) で \(V_{i}\) に属する頂点から \(V_{j}\) に属する頂点への最短路の長さの最小値を表すことにします: \[ \delta(i, j) = \min \lbrace \mathit{dist}(v_{i}, v_{j})\ |\ v_{i} \in V_{i} \text{ かつ } v_{j} \in V_{j} \rbrace \] 全ての \(i, j\) について \(\delta(i, j)\) を計算するアルゴリズムを説明してください。満点のためには実行時間が \(O(VE + kV\log V)\) であることが求められます。

  5. この問題はあなたが、そうですあなたが、ウォールストリートで職を得て大規模な経済崩壊を起こすことができるかどうかを試します! アービトラージビジネスとは、為替交換の鞘を利用して儲ける手法を言います。例えば 1 US ドルが 120 円で、1 円が 0.01 ユーロで、1 ユーロが 1.2 US ドルのとき、1 ドルを持っているトレーダーはそのドルを円に、円をユーロに、ユーロをドルに交換することで 1.44 ドルを入手できます! この例における $ → ¥ → € → $ のような通貨のサイクルのことをアービトラージサイクル (arbitrage cycle) と言います。もちろん、アービトラージサイクルを見つけて儲けるには極限まで高速なアルゴリズムが必要です。

    \(n\) 個の異なる通貨が取引される為替市場で、その交換レートが行列 \(R[1..n, 1..n]\) で与えられるとします。任意の添え字 \(i, j\) について一単位の通貨 \(i\) が \(R[i, j]\) 単位の通貨 \(j\) と交換可能なことを表すとします (\(R[i, j] \cdot R[j, i] = 1\) が成り立つとは限りません)。

    1. アービトラージサイクルが存在しないと仮定します。一単位の通貨 \(1\) から交換をしていったときに最大で入手できる通貨 \(i\) の金額が \(V[i]\) であるような配列 \(V[1..n]\) を計算するアルゴリズムを説明してください。

    2. 与えられた為替交換レートを表す行列にアービトラージサイクルが含まれるかどうかを判定するアルゴリズムを説明してください。

    3. (b) のアルゴリズムを改変してアービトラージサイクルが存在する場合にはそれを返すようにしてください。

  6. Morty は クラックスパイア迷宮 (Clackspire Labyrinth) から安定化したプランバス (plumbus) を取ってこようとしています。Morty は Rick の異次元ポータル銃を使って迷宮に入り、迷宮をプランバスの所まで進み、さらにフリーブ (fleeb) の所まで進んでプランバスを安定化させ、そして入ってきたポータルに戻って来ます。プランバスはフリーブのジュースで安定化させることができ、このジュースはフリーブを巣穴から動かすとすぐに放出されます。また安定化されていないプランバスを保管場所から 137 flink 以上遠くに動かすと爆発します。クラックスパイア迷宮はおならのようなにおいが充満しているので、Morty はなるべく短い時間で仕事を終わらせようと思っています。

    Morty は Rick からクラックスパイア迷宮の詳細な地図をもらいました。この地図は辺に非負の重みが付いた有向グラフ \(G = (V, E)\) であり、重みが距離 (単位は flink) を表します。また \(P \subset V\) と \(F \subset V\) も与えられ、\(P\) がプランバスの保管場所を、\(F\) がフリーブの巣穴を表します。Morty は探索を始める頂点 \(s \in V\) と使用するプランバスの保管場所 \(p \in P\) とフリーブの巣穴 \(f \in F\) を決めるのですが、そのとき \(p\) から \(f\) への最短路の長さは 137 flink 以下で、かつ \(s \rightsquigarrow p \rightsquigarrow f \rightsquigarrow s\) の長さは最短である必要があります。

    Morty の問題を (ゲフッ) 解くアルゴリズ (グフッ) ムを説明、解析してください。Morty がこの仕事を行えると仮定して構いません。

  7. 辺に重みの付いた有向グラフを \(G = (V, E)\) とし、辺の重みは任意の実数とします。

    1. どの頂点の組の間の最短路の長さも変化しないようにいずれかの頂点 \(v\) を取り除くことができるでしょうか? 部分グラフ \(G^{\prime} = (V \backslash \lbrace v \rbrace, E)\) であって \(G^{\prime}\) に含まれる任意の頂点の組の間の最短路の長さが \(G\) における同じ頂点の間の最短路の長さと同じであるものを \(O(V^{2})\) 時間で構築するアルゴリズムを説明してください。

    2. \(G^{\prime}\) における全ての頂点の組に対する最短路の距離が既に求まっているとします。\(v\) から他の全ての頂点への最短路の長さと他の全ての頂点から \(v\) への最短路の長さを \(O(V^{2})\) 時間で計算するアルゴリズムを説明してください。

    3. (a) と (b) を組み合わせて、 \(O(V^{3})\) 時間で動作する本文中で説明したのと異なる全組最短路アルゴリズムを説明してください (答えのアルゴリズムは Floyd-Warshall とほとんど同じです)。

  8. 次の三種類目の行列積もまた最短路アルゴリズムに関係があります。\(n \times n\) 行列 \(A, B\) の "boolean 積" あるいは"and-or 積" \(C\) は次のように定義されます: \[ C[i,j] := \bigvee_{k} \left(A[i, k] \land B[k,j] \right) \]

    1. boolean 行列積を min-plus 行列積に帰着させてください。つまり、二つの \(n \times n\) 行列の min-plus 積を \(T(n)\) 時間で計算するサブルーチン \(\textsc{MinPlusMaltiply}\) が与えられたときに、このサブルーチンを使って二つの \(n \times n\) 行列の boolean 積を \(O(T(n))\) 時間で計算するアルゴリズムを説明、解析してください。

    2. boolean 行列積を通常の行列積に帰着させてください。つまり、二つの \(n \times n\) 行列の通常の積を \(T(n)\) 時間で計算するサブルーチン \(\textsc{MatrixMaltiply}\) が与えられたときに、このサブルーチンを使って二つの \(n \times n\) 行列の通常の積を \(O(T(n))\) 時間で計算するアルゴリズムを説明、解析してください。

  9. 有向グラフ \(G\) の推移閉包 (transitive closure) とは \(G\) と同じ頂点を持つグラフであって、\(u \rightarrow v\) という辺が存在するのが \(G\) において \(u\) から \(v\) に路が存在するときに限るものを言います。

    1. ある \(2 \leq \omega < 3\) に対して、二つの \(n \times n\) 行列の boolean 行列積が \(O(n^{\omega})\) 時間で行えるとします (練習問題 9.9.(b) から \(\omega \leq 2.372864\) が分かります)。\(n\) 頂点の有向グラフの推移閉包を \(O(n^{\omega} \log n)\) 時間で計算するアルゴリズムを説明してください。

    2. \(G\) を有向非巡回グラフとします。\(G\) の推移閉包を \(O(V^{\omega})\) 時間で計算するアルゴリズムを説明してください。 [ヒント: DAG に必ず行うあの処理をして、さらに分割統治をします。\(\omega \geq 2\) を使います。]

    3. \(G\) を任意の巡回グラフとします。\(G\) の推移閉包を \(O(V^{\omega})\) 時間で計算するアルゴリズムを説明してください。 [ヒント: 任意の有向グラフを DAG にするときに必ず行うあの処理をします。]

    4. 前問の逆を行いましょう。ある \(2 \leq \alpha < 3\) について、\(n\) 頂点の有向グラフを受け取って \(O(n^{\alpha})\) 時間で推移閉包を計算するサブルーチン \(\textsc{TransitiveClosure}\) が与えられたときに、\(n \times n\) の boolean 行列積を \(O(n^{\alpha})\) 時間で計算するアルゴリズムを説明、解析してください。

  10. 次に示す再帰的なアルゴリズムが全組最短路の距離を \(O(n^{3})\) 時間で正しく計算することを示してください。問題を単純にするために、\(n\) が \(2\) のべきであると仮定してください。これまでと同じように、ヘルパー関数 \(\textsc{RecAPSP}\) に渡される配列 \(D\) は参照で渡されます。 [ヒント: このアルゴリズムは一目見ただけでは Floyd-Warshall を分かりにくくしたものに過ぎませんが、キャッシュの振る舞いがはるかに優れています8。]

    procedure \(\texttt{RecursiveAPSP}\)(\(V, E, w\))

    \(n \leftarrow |V|\)

    for \(i \leftarrow 1\) to \(n\) do

    for \(j \leftarrow 1\) to \(n\) do

    if \(i = j\) then

    \(D[i, j] \leftarrow 0\)

    if \(i \rightarrow j \in E\) then

    \(D[i, j] \leftarrow w(i \rightarrow j)\)

    else

    \(D[i,j] \leftarrow \infty\)

    \(\texttt{RecAPSP}\)(\(D, n, 1, 1, 1\))

    return \(D[1..n, 1..n]\)

    procedure \(\texttt{RecAPSP}\)(\(D,n,i,j,k\))

    if \(n=1\) then

    \(D[i,j] \leftarrow \min \lbrace D[i,j],\ D[i,k] + D[j,k] \rbrace\)

    else

    \(m \leftarrow n / 2\)

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i,\ \hphantom{,+m} j,\ \hphantom{,+m} k \hphantom{,+m} \))

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i,\ \hphantom{,+m} j,\ \hphantom{,+m} k+m \))

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i,\ \hphantom{,+m} j+m,\ k \hphantom{,+m} \))

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i,\ \hphantom{,+m} j+m,\ k+m \))

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i+m,\ j,\ \hphantom{,+m} k \hphantom{,+m} \))

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i+m,\ j,\ \hphantom{,+m} k+m \))

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i+m,\ j+m,\ k \hphantom{,+m} \))

    \(\texttt{RecAPSP}\)(\(D,\ n/2,\ i+m,\ j+m,\ k+m \))

  11. \(n\) 頂点の重み無し連結無向グラフを \(G = (V, E)\) とし、\(G\) の隣接行列を \(A[1..n, 1..n]\) とします。この問題では最短路の長さを収めた行列 \(D[1..n, 1..n]\) を計算する、Seidel によって示されたアルゴリズムを考えます。このアルゴリズムは高速な行列乗算を使っており、\(O(n^{3})\) よりも高速です。

    ある定数 \(\omega \geq 2\) に対して、二つの \(n \times n\) 行列の積を \(O(n^{\omega})\) 時間で計算するサブルーチン \(\textsc{MatrixMaltiply}\) が使えるとします。

    1. \(G\) と同じ頂点を持ち、頂点の組の間に辺があるのが \(G\) においてそれらの間に長さが 2 以下の路があるときに限られるようなグラフを \(G^{2}\) とします。一度の \(\textsc{MatrixMaltiply}\) 呼び出しと \(O(n^{2})\) 時間の追加の計算を使って \(G^{2}\) の隣接行列を計算するアルゴリズムを計算してください。

    2. \(G^{2}\) が完全グラフだとします。\(G\) の最短路を収めた行列 \(D\) を \(O(n^{2})\) 時間の追加の計算で計算するアルゴリズムを説明してください。

    3. \(G^{2}\) の最短路の距離を収めた行列 \(D^{2}\) を再帰的に計算するとします。\(G\) における \(i\) から \(j\) への最短路の長さが \(2 \cdot D^{2}[i, j]\) または \(2 \cdot D^{2}[i, j] - 1\) であることを示してください。

    4. \(G^{2}\) が完全グラフでないとします。\(X = D^{2} \cdot A\) とし、\(\deg (i)\) で \(G\) における頂点 \(i\) の次数を表すことにします。頂点 \(i\) から頂点 \(j\) への最短路の長さが \(2 \cdot D^{2}[i, j]\) となるのは \(X[i, j] > D^{2}[i,j] \cdot \deg(i)\) である場合に限ることを証明してください。

    5. \(G\) の最短路の長さを収めた行列 \(D\) を \(O(n^{\omega} \log n)\) 時間で計算するアルゴリズムを説明してください。

  12. Gideon Yuval は 1976 年、次に示す min-plus 行列乗算から通常の行列乗算への帰着を提案しました。\(A, B\) を \(n \times n\) の行列とし、要素は \(0\) から \(M\) の整数だとします。このとき \(A\) と \(B\) の min-sum 積 \[ C[i,k] = \min\limits_{j}(A[i,j] + B[j,k]) \] を全ての \(i, j\) について計算します。新しい \(n \times n\) 行列 \(A^{\prime}, B^{\prime}\) を \[ \begin{aligned} A^{\prime}[i,j] & = n^{M-A[i,j]} \\ B^{\prime}[i,j] & = n^{M-B[i,j]} \end{aligned} \] と定義し、\(A^{\prime}\) と \(B^{\prime}\) の (通常の) 積を \(C^{\prime}[i,k] = \sum_{j} A^{\prime}[i, j] \cdot B^{\prime}[j, k]\) とします。

    1. 通常の整数演算 (\(+, -, \times\)) だけを使って \(A\) から \(A^{\prime}\) を計算するアルゴリズムを説明してください。

    2. 通常の整数演算 (\(+, -, \times\)) だけを使って \(C^{\prime}\) から min-sum 積 \(C^{\prime}\) を計算するアルゴリズムを説明してください9

    3. ある \(2 \leq \omega < 3\) について、\(n \times n\) 整数行列の (通常の) 積を \(O(n^{\omega})\) 回の算術演算で計算できるとします。Yuval のアルゴリズムを使って min-sum 積 \(C\) を計算するのに算術演算は何回必要になりますか?

    4. \(n \times n\) の整数行列 \(A\) が与えられたときに、\(A\) の \(n\) "おかしな"乗を Yuval のアルゴリズムを使って求めるのに算術演算は何回必要になりますか? (\(A\) が重み付き隣接行列とすれば、その \(n\) "おかしな" 乗は最短路の長さを収めた行列であることを思い出してください)

    5. 任意の重みを持つグラフを考えたときには、Yuval のアルゴリズムが Floyd-Warshall のアルゴリズムよりも速い全組最短路アルゴリズムとはなりません。どうしてですか? ごまかしが起こっているのはどこですか?


  1. 正確には歩道 (walk) です。[return]

  2. 都市の間の距離を調べるのに印刷された道路地図を使っていた時代には、地図に三角形の "距離表" が付いていました。例えばシャンペーンとコロンバスの間の距離が知りたいときには、この表の "シャンペーン" の行と "コロンバス" の列の交点を調べます。[return]

  3. ここでも、Dijkstra のアルゴリズムの実装で使う二分ヒープをソートされていない配列に入れ替えると、全体の実行時間は (グラフに含まれる辺の数に関係なく) \(O(V^{4})\) となります。また二分ヒープをフィボナッチヒープに入れ替えると、実行時間は \(O(V(E + V\log V)) = \) \(O(VE + V^{2} \log V) = \) \(O(V^{3})\) まで落ちます。[return]

  4. ...標準的な二分ヒープの実装を使った場合は。一つ前の脚注を参照してください。[return]

  5. Shimbel は入力として距離を表す \(V \times V\) 行列を仮定していたので、彼のオリジナルのアルゴリズムの実行時間はグラフの持つ辺の数に関わらず \(O(V^{4})\) でした。[return]

  6. Raimund Seidel. On the all-pairs-shortest-path problem in unweighted undirected graphs. Journal of Computer and System Sciences, 51(3):400-403, 1995. この論文 (の少なくとも 1992 年の会議で使われたバージョン) はアルゴリズムの完全な説明と解析が論文の概要に収まっているという稀な論文の一つです。 Noga Alon, Zvi Galil, Oded Margalit*. On the exponent of the all pairs shortest path problem. Journal of Computer and System Sciences 54(2):255–262, 1997. も参照してください。[return]

  7. 発音は「cley nee (クレイニ)」です。「clean (クリーン)」、「clean-ee (クリー二)」、「clay-nuh (クレイナ)」、「dimaggio (ディマジオ)」は全て間違っています [訳注: 日本語では「クリーネ」と呼ばれますが、これも間違いです]。彼は任意の有限オートマトンに対応する正規表現が存在することを帰納法を使って証明しており、その帰納法のパターンが Floyd-Warshall の再帰方程式と本質的に同じです。[return]

  8. Joon-Sang Park, Michael Penner, and Viktor K. Prasanna. Optimizing graph algorithms for improved cache performance. IEEE Trans. Parallel and Distributed Systems 15(9):769–782, 2004. さらに広い動的計画法アルゴリズムのクラスの大幅な一般化については次の文献を参照してください: Rezaul Alam Chowdhury and Vijaya Ramachandran. Cache-oblivious dynamic programming. Proc. 17th SODA 591–600, 2006[return]

  9. 対数や割り算、それから床関数 \(\lfloor x \rfloor\) を使わないでください。信じてください ――これは開けてはいけないドアです。[return]

広告