Zenn
🏝️

閉半環を使ってグラフ上の最短距離を計算する！

2020/12/21に公開
2件
Haskell
math
math
tech

この記事は Haskell Advent Calendar 2020 21日目の記事です。

以前の記事でトロピカル行列を使ったグラフの最短経路の求め方を解説しました。

https://qiita.com/lotz/items/094bffd77b24e37bf20e

ここではトロピカルな隣接行列の累乗を収束するまで繰り返すという方法で最短経路を計算しましたが、実は閉半環という代数を考えると直接的に最短経路を求める計算が可能になります。そこで今回はその方法について解説したいと思います。

以前はHaskellのリスト [a] をベクトルとして行列を実装しましたが、今回はそれだと実装が少し煩雑になるので型レベル自然数を型引数に持つ Vector n a を中心に実装していきたいと思います。この話は以下の Functional Pearl が元になっていますが、この論文もリストを使って実装されているので Vector n a を使ってどのように実装できるかはこの記事で新しく試したところです。

http://stedolan.net/research/semirings.pdf

トロピカル半環

まずは半環とトロピカル半環を実装しましょう。前回のコメントでHaskellの Double には Infinity があるので直和型での場合分けは必要ないのではと教えてもらったので今回はそれを反映して実装してみます。

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

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

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

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

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

-- 無限大
inf :: Double
inf = 1.0 / 0.0

instance Semiring Tropical where
  oplus = min
  zero  = T inf
  otimes (T x) (T y) = T (x + y)
  one   = T 0

行列半環

半環の元を要素に持つ行列を考えると要素ごとの足し算と行列積を考えることで、行列全体が再び半環の構造を持ちます。

まずはベクトルと行列に関する基本的な演算を実装しましょう。ベクトルの実装にはvector-sizedを使います。

{-# LANGUAGE BlockArguments #-}

import GHC.TypeNats

import Data.Vector.Sized (Vector)
import qualified Data.Vector.Sized as V

-- m × n 行列
type Matrix m n a = Vector m (Vector n a)

-- | 単位行列
ident :: (KnownNat n, Semiring a) => Matrix n n a
ident = V.generate \i ->
          V.generate \j ->
            if i == j then one else zero

-- | 各成分ごとの演算
elementWise :: (a -> b -> c) -> Matrix m n a -> Matrix m n b -> Matrix m n c
elementWise op = V.zipWith (V.zipWith op)

-- | 行列の和
(!+!) :: Semiring a => Matrix m n a -> Matrix m n a -> Matrix m n a
(!+!) = elementWise oplus

-- | 内積
dot :: Semiring a => Vector n a -> Vector n a -> a
dot xs ys = V.foldr oplus zero (V.zipWith otimes xs ys)

-- | 転置
transpose :: (KnownNat m, KnownNat n) => Matrix m n a -> Matrix n m a
transpose = V.sequence

-- | 行列積
(!*!) :: (KnownNat p, KnownNat q, KnownNat r, Semiring a)
      => Matrix p q a -> Matrix q r a -> Matrix p r a
a !*! b = flip fmap a \as -> flip fmap b' \bs -> as `dot` bs
  where b' = transpose b

実装は至ってシンプルですね。

それでは作った行列を半環のインスタンスにしてみましょう。

{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}

instance (KnownNat m, Semiring a) => Semiring (Matrix m m a) where
  oplus  = (!+!)
  zero   = V.replicate (V.replicate zero)
  otimes = (!*!)
  one    = ident

これだけです。ちなみにzeroは要素が全て0（半環における足し算の単位元）であるような行列です。

実は今回ベクトルの実装に Vector n a を採用しているのはこの半環の実装をシンプルにしたかったのが大きな理由です。リスト [a] を使ったベクトル・行列でも半環のインスタンスを定義することは可能ですが、単純なリストだと型だけからサイズに関する情報が分からないので特別な値を用意する必要があります。このことが半環の実装や行列積などの関数の実装にまで波及して、場合分けを伴う煩雑な実装につながります。Vector n a は型レベルでサイズの情報を持っているおかげで実装が非常に単純になっているのです。

このあと、ブロック行列に関する操作も必要になるのでいくつか関数を実装しておきましょう。

{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}

type BlockMatrix s t u v a = ( Matrix s u a, Matrix s v a
                             , Matrix t u a, Matrix t v a )

mjoin :: BlockMatrix s t u v a -> Matrix (s+t) (u+v) a
mjoin (a, b, c, d) = V.zipWith (V.++) (a V.++ c) (b V.++ d)

msplit :: (KnownNat m, KnownNat n) => Matrix (1+m) (1+n) a -> BlockMatrix 1 m 1 n a
msplit m = ( first,  transpose right
           , transpose left, transpose rest )
  where
    (top, bottom)  = V.splitAt m
    (first, right) = V.splitAt (transpose top)
    (left, rest)   = V.splitAt (transpose bottom)

BlockMatrix は下図の様なサイズのブロック行列を表しています。

mjoin は与えられたブロック行列を結合して一つの行列にする関数です。

msplit は与えられた行列を下図のように、

  • スカラー
  • m次元の列ベクトル
  • n次元の行ベクトル
  • m×n次元の行列

からなるブロック行列に分解する関数です。なので引数の行列の型もMatrix (1+m) (1+n) aのように 1+ がついてることに注意してください。

閉半環

半環Sが以下のような関係を満たす関数(\cdot)^{*}: S \rightarrow Sを備えている時、Sを閉半環と呼びます[1]。

Sにおける任意の元aについて(\cdot)^{*}は、

a^* = 1 + aa^* = 1 + a^*a

を満たす。

(\cdot)^{*}は閉包（closure）と呼ばれ、感覚的には

a^* = 1 + a + a^2 + a^3 + \cdots

のような無限級数を考えていると思うことができます。クリーネ閉包を知っていればこの関数が閉包と呼ばれる感覚も分かるでしょう。

例えば0より大きい範囲に限定したトロピカル半環における閉包は常に1を返す関数です（今回はグラフの最短経路の距離をトロピカル半環として扱うのでこの範囲で十分でしょう）。これはトロピカル半環における閉包の満たすべき式が、 \mathbb{R}\cup\{\infty\} の世界における記号で

0 = {\rm min}(0, a + 0)

と書けることを意味しており、aが0より大きいので常に成り立つことが分かります。

これを実装してみましょう。

-- 閉半環
class Semiring a => ClosedSemiring a where
  closure :: a -> a

instance ClosedSemiring Tropical where
  closure _ = one

グラフの最短経路

今度は行列半環に閉半環の構造について考えてみましょう。

行列Aに対して閉包は

A^*=I+A+A^2+A^3+\cdots

のように思えるのでした。ここでAがあるグラフの隣接行列だとするとトロピカル半環を要素に持つときのA^nは各頂点間をnステップで到達できる最短距離になります。+は小さい方を選択する演算なので、結局A^*は各頂点間の最短距離を表す行列になっているはずです。つまり行列の閉半環を考えれば、隣接行列の閉包を考えることで最短距離が計算できるというわけです。

少し天下り式ですが、閉半環の元を成分に持つ行列Mが

M = \left( \begin{matrix} A & B \\ C & D \\ \end{matrix} \right)

のようなブロック行列として与えられている時、M^*を

\begin{matrix} M^* &=& \left( \begin{matrix} A^* + A^*B\Delta^*CA^* & A^*B\Delta^* \\ \Delta^*CA^* & \Delta^* \\ \end{matrix} \right) \\ \Delta &=& D + CA^*B \end{matrix}

と定義すると (\cdot)^* は閉包の条件を満たします。この定義はA, \Deltaの閉包に依存していますが、それぞれのサイズはMよりも小さくなっているので再帰的に計算することが可能です。Mが一行一列の行列（つまりスカラー）であれば単純に要素の閉包を取ることとします。

実際

\begin{matrix} I + MM^* &=& \left( \begin{matrix} I & O \\ O & I \\ \end{matrix} \right) + \left( \begin{matrix} A & B \\ C & D \\ \end{matrix} \right) \left( \begin{matrix} A^* + A^*B\Delta^*CA^* & A^*B\Delta^* \\ \Delta^*CA^* & \Delta^* \\ \end{matrix} \right) \\ &=& \left( \begin{matrix} I & O \\ O & I \\ \end{matrix} \right) + \left( \begin{matrix} AA^* + AA^*B\Delta^*CA^* + B\Delta^*CA^* & AA^*B\Delta^* + B\Delta^* \\ CA^* + CA^*B\Delta^*CA^* + D\Delta^*CA^* & CA^*B\Delta^* + D\Delta^* \\ \end{matrix} \right) \\ &=& \left( \begin{matrix} I + AA^* + AA^*B\Delta^*CA^* + B\Delta^*CA^* & AA^*B\Delta^* + B\Delta^* \\ CA^* + CA^*B\Delta^*CA^* + D\Delta^*CA^* & I + CA^*B\Delta^* + D\Delta^* \\ \end{matrix} \right) \\ &=& \left( \begin{matrix} (I + AA^*) + (I + AA^*)B\Delta^*CA^* & (I + AA^*)B\Delta^* \\ CA^* + (D + CA^*B)\Delta^*CA^* & I + (D + CA^*B)\Delta^* \\ \end{matrix} \right) \\ &=& \left( \begin{matrix} A^* + A^*B\Delta^*CA^* & A^*B\Delta^* \\ CA^* + \Delta\Delta^*CA^* & I + \Delta\Delta^* \\ \end{matrix} \right) \\ &=& \left( \begin{matrix} A^* + A^*B\Delta^*CA^* & A^*B\Delta^* \\ \Delta^*CA^* & \Delta^* \\ \end{matrix} \right) \\ &=& M^* \end{matrix}

となり閉包の条件を満たしていることが分かります（I+M^*Mも同様）。

さっそく行列の閉包をHaskellで実装してみましょう。と言いたいところなのですが、実はGHCで定義されている型レベル自然数Natは再帰的な定義になっていないため[2]、定義通りに実装することが困難です。このことについて Haskell-jp Slack の #question チャネルで質問したところ @mod_poppo さんから明快な回答を頂き問題を回避した実装を教えていただけたので今回はそちらを使った実装を載せたいと思います（Haskell-jp Slackへの登録はこちらから）。

ポイントは closure :: Matrix m m a -> Matrix m m a を実装する時に、mは0であるか？そうでなければ何らかの型レベル自然数nが存在してm ~ 1+nとなるか（そうでないとmsplitが使えない）ということを型推論させる必要があるということです。そのために以下のようなデータ型と関数を用意します。

{-# LANGUAGE GADTs #-}
{-# LANGUAGE KindSignatures #-}

import Data.Proxy
import Data.Type.Equality
import Unsafe.Coerce (unsafeCoerce)

data NatCons (n :: Nat) where
  Zero :: NatCons 0
  Succ :: (KnownNat n, m ~ (1 + n)) => Proxy n -> NatCons m

{-# NOINLINE natCons #-}
natCons :: KnownNat n => Proxy n -> NatCons n
natCons proxy = case sameNat proxy (Proxy :: Proxy 0) of
                  Just Refl -> Zero
                  Nothing -> case someNatVal (natVal proxy - 1) of
                                SomeNat proxy' -> unsafeCoerce (Succ proxy')

GADTsの力によりnatConsの結果がZeroであればm ~ 0であり、Succ proxyであればm ~ 1+n となっていることが保証されます。これを使えば行列半環の閉包は以下のように実装することができます。

instance (KnownNat m, ClosedSemiring a) => ClosedSemiring (Matrix m m a) where
  closure m = case natCons (V.length' m) of
    Zero       -> m
    Succ proxy ->
      let (a, b, c, d) = msplit m
          a' = fmap (fmap closure) a
          d' = closure (d !+! (c !*! a' !*! b))
       in mjoin ( a' !+! (a' !*! b !*! d' !*! c !*! a'), a' !*! b !*! d'
                ,                      d' !*! c !*! a',               d' )

実装した閉包を使って実際に最短経路の距離が計算できるのか実験してみましょう。前回の記事と同じ以下のようなグラフを考えます。

あと計算する前に、便利関数としてリストから行列を作る関数と行列をきれいに表示する関数を作っておきましょう。

import Data.Maybe (fromJust)

fromList :: (KnownNat m, KnownNat n) => [[a]] -> Matrix m n a
fromList = fromJust . V.fromList . fmap (fromJust . V.fromList)

display :: Show a => Matrix m n a -> IO ()
display = mapM_ print . V.toList . (fmap V.toList)

それでは実験を始めましょう。まずは上図グラフの隣接行列を定義します。

> a = fromList [ [T inf, T 2,   T 4,   T inf]
|              , [T inf, T 0,   T 1,   T 9  ]
|              , [T inf, T inf, T inf, T 5  ]
|              , [T 3,   T inf, T inf, T inf]
|              ] :: Matrix 4 4 Tropical
> display a
[T Infinity,T 2.0,T 4.0,T Infinity]
[T Infinity,T 0.0,T 1.0,T 9.0]
[T Infinity,T Infinity,T Infinity,T 5.0]
[T 3.0,T Infinity,T Infinity,T Infinity]

いい感じですね。それでは閉包を計算してみましょう。

> display $ closure a
[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]

同じノードに留まるコストを０と考えていることを考慮すれば、見事最短経路の距離が計算できていることが分かります 👏 便利ですね！閉半環！

＼読んでいただきありがとうございました！／
この記事が面白かったら いいね♡ をいただけると嬉しいです☺️
バッジを贈っていただければ次の記事を書くため励みになります🙌

脚注

  1. Lehmann, Daniel J. "Algebraic structures for transitive closure." Theoretical Computer Science 4.1 (1977): 59-76. ↩︎

  2. GHCの型レベル自然数を理解する - Qiita ↩︎

Discussion

YAMAMOTO YujiYAMAMOTO Yuji

「（Haskell-jp Slackへの登録はこちらから）」のリンクが間違っているようです！
https://zenn.dev/lotz/articles/(https://haskell.jp/signin-slack.html) となぜか余分な文字列が着いちゃってます！

lotzlotz

あー何故か余計な括弧がついちゃってましたね 😅
修正しました。ありがとうございます🙏