動的計画法を実現する代数〜トロピカル演算でグラフの最短経路を計算する〜

トロピカル半環と呼ばれる代数構造上のトロピカル行列を利用すると動的計画法を使ってグラフの最短経路の距離を計算するという問題が単純な行列積で解けてしまうらしい。そんな噂¹²を聞きつけて我々はその謎を解き明かすべく南国（トロピカル）の奥地へと向かった。

トロピカルな世界に行くためにはまずは代数を知る必要がある。要するに群・環・体の話だ。しかしこの記事の目的は代数学入門ではないので詳しい話は他の記事³に譲るとし、さっそく半環という概念を導入する。それは

半環は以下の性質を満たす二つの二項演算、即ち加法（和）"$+$" と乗法（積）"$\cdot$" とを備えた集合$R$を言う

1. $(R, +)$ は単位元 $0$ を持つ可換モノイドを成す:

• $(a + b) + c = a + (b + c)$
• $0 + a = a + 0 = a$
• $a + b = b + a$

1. $(R, \cdot)$ は単位元$1$を持つモノイドを成す:

• $(a \cdot b)\cdot c = a \cdot (b \cdot c)$
• $1 \cdot a = a \cdot 1 = a$

1. 乗法は加法の上に分配的である:

• $a \cdot (b + c) = (a \cdot b) + (a \cdot c)$
• $(a + b) \cdot c = (a \cdot c) + (b \cdot c)$

1. $0$-倍は$R$を零化する:

• $0 \cdot a = a \cdot 0 = 0$

出典: 半環 - Wikipedia

というもの⁴。環との違いは加法の逆元の公理がないところである。例えば整数は通常の足し算と掛け算を考えると半環になる⁵。その理由は上の性質一つ一つを当てはめて確かめてみれば良い。プログラミングで言えばダックタイピング⁶みたいなもので半環の公理を満たせば何でも半環なのだ。

ところでモノイドというワードも出てきたので紹介しておくと、

集合$S$とその上の二項演算$\cdot: S \times S \rightarrow S$ が与えられ、以下の条件
【結合律】
$S$の任意の元$a, b, c$に対して、$(a \cdot b) \cdot c = a \cdot (b \cdot c)$
【単位元の存在】
$S$の元$e$が存在して、$S$の任意の元$a$に対して$e \cdot a = a \cdot e = a$
を満たすならば、組$(S, \cdot, e)$をモノイドという。

出典: モノイド - Wikipedia

という概念。整数と足し算はモノイドになるし整数と掛け算もモノイドになる。まぁ整数は半環になるって前に言ってるからそりゃそうなんだけど。

そろそろ整数以外の例⁷も出しておくと例えばリストはモノイドになる。リストの場合、二項演算は結合(append, ++)で単位元は空リスト[]。リストがそうであるようにモノイドはプログラミングのいたるところで顔を出すのだが話し出すとキリがないので参考記事を挙げるにとどめておく（Haskellerのためのモノイド完全ガイド | 雑記帳）。またプログラミングにおける代数的構造の例は @aiya000 さんの講演が詳しいし分かりやすい（Haskellの代数的構造入門　半群・モノイド・環とは何か？ - ログミーTech）。

あと大事なのは二項演算は普通の足し算と掛け算だとは限らないということ（例えばリストでは結合だった）。代数を考える時は先入観に囚われすぎずに物事を抽象的に考えることが重要。ルール（公理）を守ってさえいれば何でもモノイドだし半環である。これから向かう先はトロピカルな世界であるので先入観に囚われすぎず、細かいことは全部忘れよう。

ようこそ！トロピカルな世界へ

これまでに何度も出てきている"トロピカル"というワードは数学の一分野であるトロピカル幾何学に由来する。そもそもなんで"トロピカル"幾何学と呼ばれるかというと、どうやらこの分野を研究し始めたのがブラジルの計算機科学者だから⁸らしい。トロピカル幾何学は代数幾何学⁹や可積分系の超離散化¹⁰や代数統計学¹¹などなど幅広い応用先を持っており話題に尽きないのだが、この記事ではトロピカル幾何学の幾何学的な話はあまり出てこないので深入りはしない。主役は冒頭にも出てきたトロピカル半環である。

${\mathbb R}\cup\{\infty\}$ に以下の演算

• $x\oplus y=\min(x,y)$
• $x\otimes y=x+y$

を入れたものをトロピカル半環（あるいはmin-plus代数）と呼ぶ。
ここで

• $\infty \oplus x = x \oplus \infty = x$
• $\infty \otimes x = x \otimes \infty = \infty$

と定義する。
これにより$\infty$は零元（加法の単位元）となる。
乗法の単位元は$0$である。

ここで${\mathbb R}$は実数であり、$\infty$は実数には含まれない何かである¹²。要するに半環における"足し算"と"掛け算"として「二元の最小値を取る操作」と「普通の足し算」を採用したのがトロピカル半環というわけだ。"掛け算"が足し算になってるので本当にまぎらわしい。ちなみに"足し算"として最小値ではなく最大値を採用してトロピカル半環と呼ぶケースもあるので文献を読む時はどちらで定義されているかに注意しよう。このなんとも奇妙な"足し算"と"掛け算"を指してトロピカル演算と呼ぶのだが、本当に半環になっているのかいくつか試してみよう。

まず"足し算"、つまり「最小値を取る操作」が結合則を満たしているのか確かめたい。もしこれを満たしていなければモノイドにならない。しかしこれは

\begin{matrix}
&& (x \oplus y) \oplus z &=& x \oplus (y \oplus z) \\
&\Leftrightarrow& \min(\min(x, y), z) &=& \min(x, \min(y, z))
\end{matrix}

が成り立つことを意味しており、これは明らかに成り立つ。

あと気になるのは分配法則。これはつまり

\begin{matrix}
&& x \otimes (y \oplus z) &=& (x \otimes y) \oplus (x \otimes z) \\
&\Leftrightarrow& x + \min(y, z) &=& \min(x + y, x + z)
\end{matrix}

が成り立つことを要求しているがこれも明らかに成り立つ。

どうやらトロピカル演算は本当に半環になっているらしい（他にも気になる条件があれば直接確かめてみよ）。

よくある笑い話にカメラマンが「1足す1はー？」と唱えると数学者が「考えている代数的構造がわからないので答えられない」という笑い話があるが¹³、トロピカル半環において「1足す1はー？」と聞かれたら1である。なぜなら、

1 \oplus 1 = \min(1, 1) = 1

となるからだ。これは1だけでなく任意の元について成立し$x\oplus x = x$、この事実を指してトロピカル半環は冪等半環(dioid)であるとも言われる。

トロピカル半環に慣れるため、いくつか計算問題を解いてみよう。

\begin{matrix}
1 \otimes (2 \oplus 3) \\
(1 \oplus 2) \otimes (3 \oplus 4) \\
0 \oplus 1 \oplus \infty \\
\end{matrix}

答えはここをクリックすると表示される

答え: 上から順番に3, 4, 0

トロピカル半環は半環であるので足し算はモノイドであり一般に逆元を持たないことには注意が必要である¹⁴。他にもトロピカル多項式など書き足りない内容は山程あるが、今回は使わないのでここでは述べず先を急ぐ事にする。

トロピカル線形代数

ここからは後にアルゴリズムを考える時に必要となるトロピカル行列を導入する。しばらく定義が続くのでまるでジャングルを進むような感覚を覚えるかもしれないが、長居はしないので安心して付いてきてほしい。

トロピカル半環の元を横に並べるとトロピカルなベクトルが得られる¹⁵。

(x, y)

例として二元並べたものを用いるがもちろん何元並べても構わない。このベクトルにトロピカル半環の元を掛ける操作を成分ごとの掛け算、ベクトルの足し算を成分同士の足し算として定める。

\begin{matrix}
\alpha \otimes (x, y) &=& (\alpha\otimes x, \alpha\otimes y) \\
(a, b) \oplus (x, y) &=& (a\oplus x, b\oplus y)
\end{matrix}

ベクトル同士の内積もトロピカル演算を使って定めておく。

(a, b) \cdot (x, y) = a \otimes x \oplus b \otimes y

トロピカル半環の元を縦横に並べればトロピカル行列を作ることが出来る。

\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22} \\
\end{pmatrix}

トロピカル半環の元との掛け算、トロピカル行列同士の足し算はベクトルの時と同様に定義する。さらにベクトルの内積を利用してトロピカル行列同士の積を定める。

\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22} \\
\end{pmatrix}

\otimes

\begin{pmatrix}
b_{11} & b_{12} \\
b_{21} & b_{22} \\
\end{pmatrix}

=

\begin{pmatrix}
a_{11} \otimes b_{11} \oplus a_{12} \otimes b_{21} & a_{11} \otimes b_{12} \oplus a_{12} \otimes b_{22} \\
a_{21} \otimes b_{11} \oplus a_{22} \otimes b_{21} & a_{21} \otimes b_{12} \oplus a_{22} \otimes b_{22} \\
\end{pmatrix}

この行列積にはdistance productという名前もついていて¹⁶、後に述べるアルゴリズムにおいて中心的な役割を果たすことになる。

トロピカル行列における単位行列を考えてみよう。通常の単位行列は

\begin{pmatrix}
1 & 0 \\
0 & 1 \\
\end{pmatrix}

のように対角線上に乗法の単位元を並べ、非対角成分には零元を並べたものであった。同様に対角線上にトロピカル半環の乗法の単位元を並べ、非対角成分にはトロピカル半環の零元を並べた行列を考えると

\begin{pmatrix}
0 & \infty \\
\infty & 0 \\
\end{pmatrix}

となる。これが単位行列だと言われるとなんとも不思議な気持ちになるが、実際に任意の行列との積を考えてみると

\begin{matrix}

&&
\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22} \\
\end{pmatrix}

\otimes

\begin{pmatrix}
0 & \infty \\
\infty & 0 \\
\end{pmatrix}

\\

&=&

\begin{pmatrix}
a_{11} \otimes 0 \oplus a_{12} \otimes \infty & a_{11} \otimes \infty \oplus a_{12} \otimes 0 \\
a_{21} \otimes 0 \oplus a_{22} \otimes \infty & a_{21} \otimes \infty \oplus a_{22} \otimes 0 \\
\end{pmatrix}

\\

&=&

\begin{pmatrix}
a_{11} \oplus \infty & \infty \oplus a_{12} \\
a_{21} \oplus \infty & \infty \oplus a_{22} \\
\end{pmatrix}

\\

&=&

\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22} \\
\end{pmatrix}

\end{matrix}

となってちゃんと単位行列の役割を果たしていることが分かる。

単位行列はもちろんトロピカル行列の行列積における単位元である。また成分が全て$\infty$であるような行列を考えるとこれはトロピカル行列の加法における単位元となる。実はトロピカル行列はこれらの加法と乗法のもとに半環になるのである。分配法則と零化の公理を満たしているかどうかは簡単なので手を動かして確かめてみると良い。

最後にトロピカル行列の行列式を紹介する。通常の行列式は対称群${\mathfrak S}_n$¹⁷を用いて

\sum_{\sigma \in {\mathfrak S}_n} {\rm sgn}(\sigma)\prod_{i=1}^n a_{i\ \sigma(i)}

と定義されるが、トロピカル半環は逆元を持つとは限らないので${\rm sgn}(\sigma)$が現れるのは都合が悪い（${\rm sgn}(\sigma)$は置換の符号であり$\pm 1$の値を取る）。なのでトロピカル行列の行列式として通常の行列ではパーマネントと呼ばれる値を採用する。

\sum_{\sigma \in {\mathfrak S}_n} \prod_{i=1}^n a_{i\ \sigma(i)}

要するに都合の悪い符号を単純に消してしまっただけである。ついでに記号もトロピカル演算のものに変えておこう。

\bigoplus_{\sigma \in {\mathfrak S}_n} \bigotimes_{i=1}^n a_{i\ \sigma(i)}

2行2列の場合に具体的に計算してみよう。

{\rm det}
\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22} \\
\end{pmatrix}

=

\bigoplus_{\sigma \in {\mathfrak S}_2} \bigotimes_{i=1}^2 a_{i\ \sigma(i)}
=

a_{11} \otimes a_{22} \oplus a_{12} \otimes a_{21}

これはつまり行と列がかぶらないように要素を取り出して合計したものの中で一番小さい値を表している。このような解を求める問題は割当問題と呼ばれている。

行列式を求めるアルゴリズムはガウス掃き出し法を使えば$O(n^3)$で計算できるが、通常のパーマネントを計算するアルゴリズムは#P完全であり多項式時間で計算することは絶望的とされている。しかしトロピカル行列の行列式は割当問題とみなすことができるのでハンガリアン法¹⁸を使えば$O(n^3)$で計算することが可能である¹⁹。

トロピカル行列の逆行列や固有値固有ベクトルなど、まだまだ道は続いているので興味がある人は以下の文献を読んでみることをオススメする。

• 経営学における代数学 -線形代数からMax-Plus代数へ-
• （上級者向け）『トロピカル幾何学入門』

トロピカルなグラフ理論

さて、いよいよグラフの最短経路の距離を行列積で計算してみよう。以下ではグラフとして重み付き有向グラフを考える。これはグラフなのでループがあっても良い。あと重みは負を許してもほとんど問題ないが細かい注釈をつけるのが面倒なのでここでは非負としておく。例えば以下のようなグラフを考える。

重み付き有向グラフの頂点iから頂点jに向かって重みwの辺が存在する場合にi行j列をwとし、辺が存在しない場合は値が∞であるような行列を重み付き有向グラフの隣接行列と呼ぶ。例えば上のグラフの隣接行列は、

\begin{pmatrix}
\infty & 2 & 4 & \infty \\
\infty & 0 & 1 & 9 \\
\infty & \infty & \infty & 5 \\
3 & \infty & \infty & \infty \\
\end{pmatrix}

となる。各行各列は頂点a,b,c,dの順に対応させた。

さてこの隣接行列には∞という値が現れており如何にもトロピカル半環と相性が良さそうである。実際、次の定理が成り立つ。

【定理】

重み付き有向グラフの隣接行列$A$をトロピカル半環の元を要素とするトロピカル行列とみなす。
この時、$A$をトロピカル行列の行列積で$m$乗した行列$A^{\otimes m}$のi行j列の値は、頂点iから頂点jまでm回で到達できる経路のうち最小の距離（重みの和）を持つものの距離と対応する。

証明する前に具体例を見てみよう。先程の例における隣接行列の累乗を考えると

\begin{matrix}
A^{\otimes 2} &=&
\begin{pmatrix}
\infty & 2 & 3 & 9 \\
12 & 0 & 1 & 6 \\
8 & \infty & \infty & \infty \\
\infty & 5 & 7 & \infty \\
\end{pmatrix} \\
A^{\otimes 3} &=&
\begin{pmatrix}
12 & 2 & 3 & 8 \\
9 & 0 & 1 & 6 \\
\infty & 10 & 12 & \infty \\
\infty & 5 & 6 & 12 \\
\end{pmatrix} \\
\end{matrix}

のようになる。実際に$m$乗の値が$m$回で到達できる最短経路の距離に対応してることを確かめてほしい。$m$回以下ではなくちょうど$m$回なので一度実数値となった成分も累乗を進めれば再び$\infty$になっている場合があることも特筆すべき点である。

それでは定理を証明してみよう。

【証明】

$m=1$の時は隣接行列そのものであるが、これが1回でいける経路の最短距離を表していることは明らか。
$m=r$の時に$A^{\otimes r}$のi,j成分を$a^{(r)}_{ij}$と置きこれがr回で到達できる経路の最短距離を表しているとすると、$m=r+1$の時の$A^{\otimes r+1}$における$a^{(r+1)}_{ij}$は行列積の定義より

a^{(r+1)}_{ij} = \bigoplus_k a^{(r)}_{ik}\otimes a_{kj} = \min\{a^{(r)}_{ik} + a_{kj} : k=1,\dots,n\}

となるが、これはiからjまでちょうどr+1回でたどり着ける経路の最短距離に他ならない。よって帰納法より示された □

トロピカル演算の和を取る操作、すなわち最小値を取る操作は、このアルゴリズムが動的計画法として効率よいものであるためにとても重要な役割を果たしているのである。

レッツ、トロピカルコーディング！

それではアルゴリズムが手に入ったので実際にコーディングしてみよう。まずはトロピカル半環を実装していく。

-- | トロピカル半環
data Tropical = T Double
              | Infty
              deriving (Show, Eq, Ord)

トロピカル半環は集合として${\mathbb R}\cup\{\infty\}$であり、実装はそっくりそのままである²⁰。次に半環を定義しよう。

-- | 半環
class Semiring a where
  -- | 加法
  oplus :: a -> a -> a

  -- | 加法の単位元
  zero :: a

  -- | 乗法
  otimes :: a -> a -> a

  -- | 乗法の単位元
  one :: a

instance Semiring Tropical where
  oplus = min

  zero = Infty

  otimes _ Infty = Infty
  otimes Infty _ = Infty
  otimes (T x) (T y) = T (x + y)

  one = T 0

半環は型クラスとして定義した。もちろんTropicalはSemiringのインスタンスとなる。次にトロピカル線形代数を実装していこう。

-- | ベクトル
type Vector a = [a]

-- | ベクトルの次元
dim :: Vector a -> Int
dim = length

ベクトルは単純にリストとして実装する。

-- | 行列
type Matrix a = [Vector a]

-- | n次元の単位行列
ident :: Semiring a => Int -> Matrix a
ident n = reverse
        . take n . map (take n)
        . tails . cycle
        $ replicate (n-1) zero ++ [one]

-- | 行列の和
matrixPlus :: Semiring a => Matrix a -> Matrix a -> Matrix a
matrixPlus = zipWith (zipWith oplus)

もし行列の型にサイズの情報まで持たせるような実装にすれば、前述の通りトロピカル行列は半環であったのでSemiring a => Matrix aをSemiringのインスタンスにすることもできる²¹²²。しかし今回はそこまでしてもこの後の実装は大して楽にならないので、簡単に行列をベクトルのリストとして実装する。

-- | 内積
dot :: Semiring a => Vector a -> Vector a -> a
dot xs ys = foldr oplus zero $ zipWith otimes xs ys

-- | 行列積
distanceProduct :: Semiring a => Matrix a -> Matrix a -> Matrix a
distanceProduct a b = [[as `dot` bs | bs <- transpose b] | as <- a]

-- | 行列の累乗
power :: Semiring a => Matrix a -> Int -> Matrix a
power a 0 = ident $ dim (head a)
power a n = a `distanceProduct` power a (n-1)

ここまで定義ができれば準備は万端だ。先程の例を実際に実行してみよう。

> a = [ [Infty, T 2,   T 4,   Infty]
|     , [Infty, T 0,   T 1,   T 9  ]
|     , [Infty, Infty, Infty, T 5  ]
|     , [T 3,   Infty, Infty, Infty]
|     ]

> mapM_ print $ a `distanceProduct` a
[Infty,T 2.0,T 3.0,T 9.0]
[T 12.0,T 0.0,T 1.0,T 6.0]
[T 8.0,Infty,Infty,Infty]
[Infty,T 5.0,T 7.0,Infty]

期待通りの結果が帰ってきた。累乗も実行してみよう。

> mapM_ print $ a `power` 3
[T 12.0,T 2.0,T 3.0,T 8.0]
[T 9.0,T 0.0,T 1.0,T 6.0]
[Infty,T 10.0,T 12.0,Infty]
[Infty,T 5.0,T 6.0,T 12.0]

> mapM_ print $ a `power` 100
[T 11.0,T 2.0,T 3.0,T 8.0]
[T 9.0,T 0.0,T 1.0,T 6.0]
[T 19.0,T 10.0,T 11.0,T 16.0]
[T 14.0,T 5.0,T 6.0,T 11.0]

100乗の計算なんてのも朝飯前だ。ところで6乗の計算結果は

A^{\otimes 6} =
\begin{pmatrix}
11 & 2 & 3 & 8 \\
9 & 0 & 1 & 6 \\
19 & 10 & 11 & 16 \\
14 & 5 & 6 & 11 \\
\end{pmatrix} \\

であり100乗の計算結果と一致してるのが見て分かる。実はこの例のグラフでは隣接行列の6乗以降は全て同じ値になる。これは頂点bにとどまるコストが0であるので、ある程度回数が多くなると途中bにとどまることでコストを調節できるようになり最短経路が変わらなくなってしまうからだ。この考え方を応用して$m$回以内に到達できる最短経路の距離を計算することが出来る。それは隣接行列に単位行列を足し算して累乗を計算すれば良い。そうすれば各頂点にコスト0のループを付け加えたのと同じ計算になるので最短経路がちょうどm回でなくても頂点にとどまって調整するという選択ができるようになる。

> mapM_ print $ (a `matrixPlus` ident 4) `power` 1
[T 0.0,T 2.0,T 4.0,Infty]
[Infty,T 0.0,T 1.0,T 9.0]
[Infty,Infty,T 0.0,T 5.0]
[T 3.0,Infty,Infty,T 0.0]

> mapM_ print $ (a `matrixPlus` ident 4) `power` 2
[T 0.0,T 2.0,T 3.0,T 9.0]
[T 12.0,T 0.0,T 1.0,T 6.0]
[T 8.0,Infty,T 0.0,T 5.0]
[T 3.0,T 5.0,T 7.0,T 0.0]

> mapM_ print $ (a `matrixPlus` ident 4) `power` 3
[T 0.0,T 2.0,T 3.0,T 8.0]
[T 9.0,T 0.0,T 1.0,T 6.0]
[T 8.0,T 10.0,T 0.0,T 5.0]
[T 3.0,T 5.0,T 6.0,T 0.0]

> mapM_ print $ (a `matrixPlus` ident 4) `power` 4
[T 0.0,T 2.0,T 3.0,T 8.0]
[T 9.0,T 0.0,T 1.0,T 6.0]
[T 8.0,T 10.0,T 0.0,T 5.0]
[T 3.0,T 5.0,T 6.0,T 0.0]

今度は$A^{\otimes 3}$以降が同じ値になっていることが分かる。頂点にとどまるという選択肢が増えた場合、各頂点を1回以上通る経路は最短になりえないので必ず 頂点の数-1回 以内の経路が最短になる。なので累乗もそれだけ行えば十分なのである。実は閉半環という代数構造を利用すれば反復計算を行って収束するのを待たなくても直接最短経路を計算することもできる（『閉半環を使ってグラフ上の最短距離を計算する！』）。

以上のように特定の頂点間の最短経路を求めるのではなく全頂点間の最短経路を求める問題は全点対間最短経路問題と呼ばれており、トロピカル行列の行列積による解法は実はワーシャル-フロイド法というよく知られたアルゴリズムに一致している。²³

経路にもトロピカルな代数を実装しよう

これまでのアルゴリズムでは最短経路の距離しかわからなかった。しかしどうせなら最短経路が何であるかをちゃんと計算したい。そこで経路自体にトロピカル演算を考えることで最短経路を同じアルゴリズムを使って計算してみる。

まずは経路を扱うためにグラフをちゃんと実装する。

-- | グラフの頂点
type Node = Char

-- | グラフの辺
type Edge = (Node, Node)

-- | 重み
type Weight = Tropical

-- | グラフの辺に重みを対応付ける関数
weight :: Edge -> Weight
weight ('a', 'b') = T 2
weight ('a', 'c') = T 4
weight ('b', 'b') = T 0
weight ('b', 'c') = T 1
weight ('b', 'd') = T 9
weight ('c', 'd') = T 5
weight ('d', 'a') = T 3

今回は重みも付いているので辺と重みを対応付ける関数も用意している。

-- | グラフ上の経路
data Path = P [Edge]
          | Prohibited
          deriving Show

-- | 経路の重みを計算する
pathWeight :: Path -> Weight
pathWeight Prohibited = Infty
pathWeight (P es)     = foldr otimes one $ map weight es

経路は辺のリストと禁止されている経路を表すProhibitedの直和として定義する。経路の重みは経路に含まれる辺の重みを全て足したものとする。この禁止路まで含めた経路には半環の構造を入れることが出来る。

minOn :: Ord b => (a -> b) -> a -> a -> a
minOn f a b
  | f a <= f b = a
  | otherwise  = b

instance Semiring Path where
    oplus = minOn pathWeight

    zero = Prohibited

    otimes _ Prohibited  = Prohibited
    otimes Prohibited _  = Prohibited
    otimes (P xs) (P ys) = P $ xs ++ ys

    one = P []

実はpathWeightは経路の半環からトロピカル半環への準同型写像になっている。

さて最短経路を求めるために隣接行列ではなく1ステップの経路を並べた次のような行列を定義する。

-- | 経路行列
pathMatrix :: Matrix Path
pathMatrix = [ [Prohibited,     P [('a', 'b')], P [('a', 'c')], Prohibited]
             , [Prohibited,     P [('b', 'b')], P [('b', 'c')], P [('b', 'd')]]
             , [Prohibited,     Prohibited,     Prohibited,     P [('c', 'd')]]
             , [P [('d', 'a')], Prohibited,     Prohibited,     Prohibited]          

隣接行列の時と同様に単位行列を足し算して累乗してみよう。行列に関する演算は多相的に作っていたので再実装する必要がないのは本当に楽である。

> mapM_ print $ (pathMatrix `matrixPlus` ident 4) `power` 3
[P [],P [('a','b')],P [('a','b'),('b','c')],P [('a','b'),('b','c'),('c','d')]]
[P [('b','c'),('c','d'),('d','a')],P [('b','b'),('b','b'),('b','b')],P [('b','b'),('b','b'),('b','c')],P [('b','b'),('b','c'),('c','d')]]
[P [('c','d'),('d','a')],P [('c','d'),('d','a'),('a','b')],P [],P [('c','d')]]
[P [('d','a')],P [('d','a'),('a','b')],P [('d','a'),('a','b'),('b','c')],P []]

おっと、このままでは同じ場所にとどまる自明な経路も含んでしまっているので、そういうものは削除するようにしよう。

-- | 自明な経路を削除する
normalize :: Path -> Path
normalize Prohibited = Prohibited
normalize (P es) = P $ filter (\e -> weight e /= one) es

> mapM_ print $ map (map normalize) $ (pathMatrix `matrixPlus` ident 4) `power` 3
[P [],P [('a','b')],P [('a','b'),('b','c')],P [('a','b'),('b','c'),('c','d')]]
[P [('b','c'),('c','d'),('d','a')],P [],P [('b','c')],P [('b','c'),('c','d')]]
[P [('c','d'),('d','a')],P [('c','d'),('d','a'),('a','b')],P [],P [('c','d')]]
[P [('d','a')],P [('d','a'),('a','b')],P [('d','a'),('a','b'),('b','c')],P []]

これで各頂点間の最短経路が求まった。距離が知りたければpathWeightをmapすれば良い。

これまでのコードの実装をRepl.itにも用意しておいたので参考にして欲しい。
https://repl.it/@lotz84/TropicalGraphPath

旅の終わりに

トロピカル半環と呼ばれる代数構造上のトロピカル行列を利用すると動的計画法を使ってグラフの最短経路の距離を計算するという問題が単純な行列積で解けてしまう。この噂は本当だった。さらにグラフ上の経路にトロピカル演算を実装することによって最短経路そのものまで求めることも出来た。トロピカル半環を使えば他にも整数計画法の問題¹⁹やNeedleman-Wunschアルゴリズムと呼ばれるDNA解析に使われる動的計画法を使ったアルゴリズムに応用できる²⁴ことも知られている。まだまだ他にもトロピカル半環だと思うことで解釈できるアルゴリズムはたくさんあるだろう。そうすることで計算量の評価が良くなったり応用の幅が広がる可能性がある。将来エンティティに備わる自然なトロピカル半環の構造の上で処理を書くと動的計画法を使った効率のいい計算が勝手に実装される世界が来るかもしれない。そんなプログラミングの未来に思いを馳せながらこの短いトロピカルの旅を終えることにする。

1. この噂は現代数学2019年6月号の「しゃべくり線形代数」という連載で耳にした ↩

2. Mohri, Mehryar. "Semiring frameworks and algorithms for shortest-distance problems." Journal of Automata, Languages and Combinatorics 7.3 (2002): 321-350. ↩

3. 例えば Swiftで代数学入門 〜 1. 数とは何か？ ↩

4. 4の条件っているんだ。環だったら$0 \cdot a = (0 + 0) \cdot a = 0 \cdot a + 0 \cdot a$で両辺から$0 \cdot a$の逆元を引いて$0 \cdot a = 0$を導出できるけど半環は逆元が存在するとは限らないから公理に入れる必要があるのか。 ↩

5. もちろん整数は"半"だけでなく普通の環になる。環であれば半環である。 ↩

6. 「もしもそれがアヒルのように歩き、アヒルのように鳴くのなら、それはアヒルに違いない」 ダック・タイピング - Wikipedia ↩

7. マスロフ和という和を考えると他にもいろんな例を作ることが出来る。マスロフ和は極限をとると最大値や最小値を与える操作になり後に述べるトロピカル演算へのうまい橋渡しになる。この極限を取る操作はトロピカル化（Maslovの脱量子化）と呼ばれる。参考: マスロフ式算数がやたらに面白いんですけど - 檜山正幸のキマイラ飼育記 (はてなBlog) ↩

8. と英語のWikipediaに書いてあった Tropical geometry#History - Wikipedia ↩

9. 石川 剛郎 "トロピカル幾何学入門" ↩

10. 前野 俊昭 "トロピカルな幾何" ↩

11. Lior Pachter and Bernd Sturmfels "Tropical geometry of statistical models" ↩

12. $\infty$は任意の実数より大きな数に見えるかもしれないが（もちろんそれを意図して$\infty$という記号で書いてる）、慣れないうちはただの記号でありそれ以上でも以下でもないと考えるほうが余計な混乱をしなくて済む。$\infty$が何であるかは$\infty$と実数の間に定義された演算だけが示しているのである。 ↩

13. 全員数学者の会合で、ホテルの方に集合写真をお願いしたら「1+1は～」と聞かれた→「どうすれば数学者に2と言わせられるのか？」いろいろな策略が集まる - Togetter ↩

14. どうしても逆元が欲しい場合は対称化という操作を行うことで逆元を持つようにトロピカル半環を拡張することも出来る。 参考: マックス代数によるシステム理論の基礎 ↩

15. もっと抽象的にトロピカル半環上の加群から始めて基底を取って数ベクトルを作ったり、線形写像からトロピカル行列を作ったりはできるのだろうか？ ↩

16. Min-plus matrix multiplication - Wikipedia ↩

17. 対称群は置換全体のなす群: 対称群 - Wikipedia ↩

18. 割当問題のハンガリアン法をpythonで実装してみた - Qiita ↩

19. D. Maclagan, B. Sturmfels "Introduction to Tropical Geometry" ↩ ↩²

20. コメントで指摘していただいたようにHaskellのDoubleには$\infty$に対応する概念が含まれているので、それを利用すればDoubleをトロピカル半環として直接Semiringのインスタンスにすることもできる。 ↩

21. https://twitter.com/lotz84_/status/1147702560792313856 ↩

22. 対角成分に同じ要素を持つ行列の値を特別に作ることで型レベルに行列のサイズを持たなくてもSemiringのインスタンスを作ることができる。Fun with Semirings ↩

23. 半環に加えてclosureと呼ばれる演算を考えることで、停止するまで累乗を繰り返すという操作を経ずとも結果が求まるアルゴリズムを考えることもできる。Fun with Semirings ↩

24. Lior Pachter and Bernd Sturmfels "Algebraic Statistics for Computational Biology" ↩