Haskell中的Floyd-Warshall算法

Question

我在研究弗洛伊德-沃肖尔算法。现在我已经成功地在Haskell中实现了它，我实现它的方式类似于在命令式语言中实现它的方式(也就是说，使用列表的列表来模拟2D数组)，但是考虑到访问列表中的元素比访问数组中的元素慢得多，这真的是效率低下。

在Haskell中有没有更聪明的方法来做到这一点？我认为我可以通过连接一些列表来做到这一点，但总是失败。

我的代码：

floydwarshall :: [[Weight]] -> [[Weight]]
floydwarshall lst = fwAlg 1 $ initMatrix 0 $ list2matrix lst

fwAlg :: Int -> [[Weight]] -> [[Weight]]
fwAlg k m | k < rows m = let n = rows m
                             m' = foldl (\m (i,j) -> updateDist i j k m) m [(i,j) | i <- [0..n-1], j <- [0..n-1]]
                        in fwAlg (k+1) m'
          | otherwise = m

-- a special case where k is 0
initMatrix :: Int -> [[Weight]] -> [[Weight]]
initMatrix n m = if n == rows m then m else initMatrix (n+1) $ updateAtM 0.0 (n,n) m

updateDist :: Int -> Int -> Int -> [[Weight]] -> [[Weight]]
updateDist i j k m =
    let w = min (weight i j m) (weight i k m + weight k j m)
    in updateAtM w (i, j) m

weight :: Vertice -> Vertice -> [[Weight]] -> Weight
weight i j m = let Just w = elemAt (i, j) m in w

Answer 1

该算法有一个规则的访问模式，因此我们可以避免大量的索引，并且仍然使用列表编写它，具有(我认为)与命令式版本相同的渐近性能。

如果您确实希望使用数组来提高速度，您可能仍然希望对行和列执行类似的批量操作，而不是读取和写入单个单元格。

-- Let's have a type for weights.  We could use Maybe but the ordering
-- behaviour is wrong - when there's no weight it should be like
-- +infinity.
data Weight = Weight Int | None deriving (Eq, Ord, Show)

addWeights :: Weight -> Weight -> Weight
addWeights (Weight x) (Weight y) = Weight (x + y)
addWeights _ _ = None

-- the main function just steps the matrix a number of times equal to
-- the node count.  Also pass along k at each step.
floydwarshall :: [[Weight]] -> [[Weight]]
floydwarshall m = snd (iterate step (0, m) !! length m)

-- step takes k and the matrix for k, returns k+1 and the matrix for
-- k+1.
step :: (Int, [[Weight]]) -> (Int, [[Weight]])
step (k, m) = (k + 1, zipWith (stepRow ktojs) istok m)
  where
    ktojs = m !! k  -- current k to each j
    istok = transpose m !! k  -- each i to current k

-- Make shortest paths from one i to all j.
-- We need the shortest paths from the current k to all j
-- and the shortest path from this i to the current k
-- and the shortest paths from this i to all j
stepRow :: [Weight] -> Weight -> [Weight] -> [Weight]
stepRow ktojs itok itojs = zipWith stepOne itojs ktojs
  where
    stepOne itoj ktoj = itoj `min` (itok `addWeights` ktoj)

-- example from wikipedia for testing
test :: [[Weight]]
test = [[Weight 0, None, Weight (-2), None],
        [Weight 4, Weight 0, Weight 3, None],
        [None, None, Weight 0, Weight 2],
        [None, Weight (-1), None, Weight 0]]

Answer 2

我不知道如何实现最佳性能，但我可以给您一些使代码抽象化的技巧，这样您就可以更容易地进行性能调优。

首先，如果您在更改数据类型时不必重写所有内容，那就太好了。现在，您已经将所有关于列表列表的内容具体化了，所以让我们看看是否可以将其抽象出来。首先，我们必须弄清楚你的最小矩阵接口是什么。看一眼您的代码，您似乎拥有initMatrix、list2matrix、rows、elemAt和updateAtM。这些是查询或修改矩阵的函数，您需要实现这些函数才能为不同的matrix类型创建此代码的新版本。

组织此接口的一种方法是创建一个类。例如：

class Matrix m where
  list2matrix :: [[a]] -> m a
  matrix2List :: m a -> [[a]]
  rows :: m a -> Int
  elemAt :: Int -> Int -> m a -> a
  updateAtM :: a -> (Int, Int) -> m a -> m a
  setDiag :: a -> m a -> m a

(我继续添加了一个用于提取结果的matrix2List函数，并将initMatrix重命名/修改为setDiag，这感觉更通用一些。)

然后，我们可以更新您的代码以使用这个新类：

floydwarshall :: Matrix m => [[Weight]] -> m Weight
floydwarshall lst = fwAlg 1 $ initMatrix $ list2matrix lst

fwAlg :: Matrix m => Int -> m Weight -> m Weight
fwAlg k m | k < rows m = let n = rows m
                             m' = foldl (\m (i,j) -> updateDist i j k m) m [(i,j) | i <- [0..n-1], j <- [0..n-1]]
                        in fwAlg (k+1) m'
          | otherwise = m

initMatrix :: Matrix m => m Weight -> m Weight
initMatrix = setDiag 0

updateDist :: Matrix m => Int -> Int -> Int -> m Weight -> m Weight
updateDist i j k m =
    let w = min (elemAt i j m) (elemAt i k m + elemAt k j m)
    in updateAtM w (i, j) m

dist :: Matrix m => Int -> Int -> Int -> m Weight -> Weight
dist i j 0 m = elemAt i j m
dist i j k m = min (dist i j (k-1) m) (dist i k (k-1) m + dist k j (k-1) m)

现在我们要做的就是开始定义一些Matrix类型，看看性能如何！

让我们从列表开始，因为您已经完成了这项工作。我们将不得不使用一个新类型的包装器来让GHC高兴，但是忽略包装和解包，这在道德上与您编写的代码是相同的：

newtype ListMatrix a = ListMatrix { getListMatrix :: [[a]] }

instance Matrix ListMatrix where
  list2matrix = ListMatrix
  matrix2List = getListMatrix
  rows = length . getListMatrix
  elemAt i j (ListMatrix m) = m !! i !! j
  updateAtM a (i,j) (ListMatrix m) =
    let (firstRows, row:laterRows) = splitAt i m
        (firstCols, _:laterCols) = splitAt j row
    in ListMatrix $ firstRows <> ((firstCols <> (a:laterCols)):laterRows)
  setDiag x = go 0
    where go n m = if n == rows m then m else go (n+1) $ updateAtM x (n,n) m

(此外，我还填写了elemAt和updateAtM。)您应该能够运行

matrix2List @ListMatrix $ floydwarshall myList

并获得与当前相同的结果(和性能)。

现在，让我们开始实验吧！我们所需要做的就是定义Matrix的新实例，看看会发生什么。也许我们应该尝试纯函数：

data FunMatrix a = FunMatrix { size :: Int, getFunMatrix :: Int -> Int -> a }

instance Matrix FunMatrix where
  list2matrix l = FunMatrix (length l) (\i j -> l !! i !! j)
  matrix2List (FunMatrix s f) = (\i -> f i <$> [0..s-1]) <$> [0..s-1]
  rows = size
  elemAt i j m = getFunMatrix m i j
  updateAtM a (i,j) (FunMatrix s f) = FunMatrix s (\i' j' -> if i==i' && j==j' then a else f i' j')
  setDiag x (FunMatrix s f) = FunMatrix s (\i j -> if i==j then x else f i j)

这是如何表现的？一个问题是，开始的查找函数仍然只是对列表列表进行索引，速度很慢。一种修复方法是先转换为数组或向量，然后再进行索引。因为我们已经很好地抽象了所有东西，所有需要更改的就是这里的list2matrix定义，您可能会得到很好的性能提升！

关于性能的主题，我可以指出另一个注意事项。dist的定义做了一些严肃的“动态编程”。如果您直接在数组中写入和读取，这可能会很好地工作，但在这种递归形式中，您可能最终会做大量重复的工作。一种解决方法是使用memoize。我的goto memoization包是MemoTrie，它使得记忆事情变得非常容易。在这种情况下，您可以将dist更改为：

dist :: Matrix m => m Weight -> Int -> Int -> Int -> Weight
dist m = go'
  where
    go' = memo3 go
    go i j 0 = elemAt i j m
    go i j k = min (go' i j (k-1)) (go' i k (k-1) + go' k j (k-1))

这可能会给你一点提振！

您可能会考虑采纳@Chi的建议并使用STUArray，但您会遇到一个问题：STUArray接口要求数组查找必须是monad。仍然可以使用我在上面展示的抽象方法，但您必须更改函数的类型。而且，由于您更改了接口中的类型，因此需要将算法代码更新为一元。这可能有点痛苦，但可能有必要获得最佳性能。