Haskellで動的計画法を実装する２つの方法

出典： Easily Solving Dynamic Programming Problems in Haskell by
Memoization of Hylomorphisms

ザ圏論的やり方としては①Dynamorphism、手続き的な方法として②STモナドが挙げられる。

DynamorphismはHylomorphismをメモ化したようなもので、詳しくはlotz氏のサイトを参照してほしい。
Haskellerとしては、Dynamorphismはとても憧れる手法である。しかし、思ったよりも速度が出ない。。
このスクラップに二通りのLCSの解法を記載したが、いずれもTLEであった。
lotz氏によると、メモ化されたデータ構造にはO(n)でしかアクセスできないことが理由とのこと。

この記事では、STモナドによるメモ化再帰を用いた動的計画問題の解法を示し、半環を用いることで多くの異なる動的計画法に対する共通性を表現できることを見ていきたい。
例としては、冒頭の出典にあるナップサック問題、編集距離、ランダムウォークなどを使用していく。

半環とは

集合 Xとその上の2つの演算⊕と⊗について、以下の条件を満たす(X,⊕,⊗)を半環(semiring) という.

(X,⊕,0)が可換モノイド

(X,⊗,1)がモノイド

分配法則
a⊗(b⊕c)=(a⊗b)⊕(a⊗c)=(b⊕c)⊗a

零化
0⊗a=a⊗0=0
出典: (半環 - Wikipedia)[https://ja.wikipedia.org/wiki/半環]

class Semiring s where
  (<+>) :: s -> s -> s
  (<.>) :: s -> s -> s
  zero :: s
  one :: s

半環の例： Max-Plus代数・　ブール代数

⊕としてmax, ⊗として+を取ると、これはトロピカル半環やMax-Plus代数と呼ばれる半環となる。
同様に⊕としてminをとっても、半環となる。また、⊕としてor, ⊗としてandをとると、ブール代数と呼ばれる半環となる。

newtype MaxPlus = MaxPlus { unMaxPlus :: Int } deriving (Eq, Ord, Show)
newtype MinPlus = MinPlus { unMinPlus :: Int } deriving (Eq, Ord, Show)

instance Semiring MaxPlus where
  t1@(MaxPlus v1) <+> t2@(MaxPlus v2) = MaxPlus (max v1 v2)
  t1@(MaxPlus v1) <.> t2@(MaxPlus v2)
    | t1 == zero = zero
    | t2 == zero = zero
    | otherwise = MaxPlus (v1 + v2)
  zero = MaxPlus minBound
  one = MaxPlus 0

instance Semiring MinPlus where
  t1@(MinPlus v1) <+> t2@(MinPlus v2) = MinPlus (min v1 v2)
  t1@(MinPlus v1) <.> t2@(MinPlus v2)
    | t1 == zero = zero
    | t2 == zero = zero
    | otherwise = MinPlus (v1 + v2)
  zero = MinPlus maxBound
  one = MinPlus 0
  
 instance Semiring Boolean where
  t1@(Boolean v1) <+> t2@(Boolean v2) = Boolean (v1 || v2)
  t1@(Boolean v1) <.> t2@(Boolean v2)
    | t1 == zero = zero
    | t2 == zero = zero
    | otherwise = Boolean (v1 && v2)
  zero = Boolean False
  one = Boolean True

動的計画法に共通する構造は半環である

動的計画法をいくつか解くと、非常によく似ている構造があることがわかる。
それは、より小さい部分問題の結果(napsackとeditDistanceの加算、randomWalkの積)とすべての結果(napsackの最大、editDistanceの最小値、randomWalkの和)の組み合わせで表現できるということだ。

具体的には、動的計画法の問題に対して、適切な半環を設定し、問題の条件を部分問題に分ける関数(subproblems)と境界での値を返す関数(isTrivial)をもつDPProblemというデータ型に変換し、下記のdpSolveを適応させる。

{-# LANGUAGE FlexibleContexts, ScopedTypeVariables #-}
import Data.Array.ST
import Control.Monad.ST

data DPProblem p sc = DPProblem {
  start :: p,
  getRange :: (p, p),
  isTrivial :: p -> Maybe sc,
  subproblems :: p -> [(sc, p)]
}

dpSolve :: forall i sc s. (Semiring sc, Ix i, Eq sc) => DPProblem i sc -> ST s sc
dpSolve dp = do
  memo <- newArray (getRange dp) zero :: ST s (STArray s i sc)
  go (start dp) memo where
    go p memo
      | Just val <- isTrivial dp p = return val
      | otherwise = do
        res <- readArray memo p
        if res /= zero then return res
        else do
          ret <- foldM (\acc (s, sp) -> (<+>) acc <$> ((<.>) s <$> go sp memo)) zero (subproblems dp p)
          writeArray memo p ret
          return ret

具体例

競プロの問題を用いて、実際の解法を示す。問題の本質は半環の選択とsubproblemsにある。

1次元DP

dungeon :: (Array Int Int,Array Int Int,  Int) -> DPProblem Int MinPlus
dungeon (as,bs, n) = DPProblem {
  start = n,
  getRange = (1,n),
  isTrivial = \p -> case p of
    1 -> Just $ MinPlus 0
    2 -> Just $ MinPlus $ as ! 2
    _ -> Nothing,
  subproblems = \p -> [(MinPlus $ as ! p, p-1), (MinPlus $ bs ! p, p-2)]
}

solve arr = runST $ dpSolve $ dungeon (as,bs,n)

部分和問題

subsetSum :: (Array Int Int, Int, Int) -> DPProblem (Int,Int) Boolean
subsetSum (as, n, k) = DPProblem {
  start = (n,k),
  getRange = ((0,0), (n,k)),
  isTrivial = \(i,j) -> 
    if  i == 0 then 
      if j == 0 then Just (Boolean True)
      else Just (Boolean False)
    else if j == 0 then Just (Boolean True)
    else if j < 0 then Just (Boolean False)
    else Nothing,
  subproblems = \(i,j) -> [(Boolean True, (i-1,j-as ! i)), (Boolean True, (i-1,j))]
}

ナップサック問題

knapsack:: (Array Int (Int,Int), Int, Int) -> DPProblem (Int,Int) MaxPlus
knapsack (wvs, n, w) = DPProblem {
  start = (w,1),
  getRange = ((0,1), (w,n)),
  isTrivial = \(c,i) -> if i > n then MaxPlus 0 then Nothing,
  subproblems = \(c,i) -> 
    let (wi, vi) = wvs ! i
    in if c >= wi then [(MaxPlus vi, (c-wi, i+1)),(MaxPlus 0,(c,i+1))] else [(MaxPlus 0, (c, i+1))]
}
solve arr = runST $ dpSolve $ knapsack arr

LCS

lcs :: (Array Int Char, Array Int Char) -> DPProblem (Int,Int) MaxPlus
lcs (s1, s2) = DPProblem {
  start = (1,1),
  getRange = ((1,1), (n1,n2)),
  isTrivial = \(i,j) -> if i > n1 || j > n2 then MaxPlus 0 then Nothing,
  subproblems = \(i,j) -> 
    let c1 = s1 ! i
        c2 = s2 ! j
    in if c1 == c2 then [(MaxPlus 1, (i+1, j+1))] else [(MaxPlus 0, (i+1, j)), (MaxPlus 0, (i, j+1))]
} where
  n1 = length s1
  n2 = length s2
  
solve arr = runST $ dpSolve $ lcs arr

編集距離

editDistance :: (Array Int Char, Array Int Char) -> DPProblem (Int,Int) MinPlus
editDistance (s1, s2) = DPProblem {
  start = (1,1),
  getRange = ((1,1), (n1+1,n2+1)),
  isTrivial = \(i,j) -> if i > n1 && j > n2 then MinPlus 0 then Nothing,
  subproblems = \(i,j) -> 
    let c1 = s1 ! i
        c2 = s2 ! j
    in 
      if i > n1 then [(MinPlus 1, (i, j+1))] 
      else if j > n2 then [(MinPlus 1, (i+1, j))]
      else if c1 == c2 then [(MinPlus 0, (i+1, j+1))] 
      else [(MinPlus 1, (i+1, j)), (MinPlus 1, (i, j+1)), (MinPlus 1, (i+1, j+1))]
} where
  n1 = length s1
  n2 = length s2

solve arr = runST $ dpSolve $ editDistance arr

区間DP

blockgame :: (Array Int (Int,Int), Int) -> DPProblem (Int,Int) MaxPlus
blockgame (pas, n) = DPProblem {
  start = (1,n),
  getRange = ((1,1), (n,n)),
  isTrivial = \(l,r) -> if l == r then MaxPlus 0 then Nothing,
  subproblems = \(l, r) ->  
    let (p1,a1) = pas ! l
        (p2,a2) = pas ! r
        v1 = if l <= p1 && p1 <= r then a1 else 0
        v2 = if l <= p2 && p2 <= r then a2 else 0
    in [(MaxPlus v1, (l+1, r)), (MaxPlus v2, (l, r-1))]
}
solve arr = runST $ dpSolve $ blockgame arr

最長部分回文(Longest Subpalindrome)

palindrome :: Array Int Char -> DPProblem (Int, Int) MaxPlus
palindrome s = DPProblem {
  start = (1,n),
  getRange = ((1,1), (n,n)),
  isTrivial = \(l,r) -> if l > r then MaxPlus 0 then Nothing,
  subproblems = \(i,j) -> 
    let c1 = s ! i
        c2 = s ! j
    in 
      if i == j then [(MaxPlus 1, (i+1, j-1))]
      else if c1 == c2 then [(MaxPlus 2, (i+1, j-1)),(MaxPlus 0, (i+1, j)), (MaxPlus 0, (i, j-1))] 
      else [(MaxPlus 0, (i+1, j)), (MaxPlus 0, (i, j-1))]
} where
  n = length s

ビットDP

-- asはクーポン券iを選んだ場合の整数表現vをarrayで表現したもの
allFree :: (Array Int Int, Int, Int) -> DPProblem (Int, Int) MinPlus
allFree (as, m,n) = DPProblem {
  start = (1,0),
  getRange = ((1,0), (m+1,bit n )),
  isTrivial = \(i,s) -> if i > m && s == bit n - 1 then MinPlus 0 then Nothing,
  subproblems = \(i,s) ->
    if i > m then [] else 
      let a = as ! i in
      if s .|. a /= s then [(MinPlus 1, (i+1, s .|. a)), (MinPlus 0, (i+1, s))] else [(MinPlus 0, (i+1, s))]
}

巡回セールスマン問題

--distsは二点間の距離、貰うDPで実装
travelingSalesman :: (Array (Int, Int) Double, Int) -> DPProblem (Int, Int) MinPlus
travelingSalesman (dists, n) = DPProblem {
  start = (bit n-1, 1),
  getRange = ((bit 0, 1), (bit n-1, n)),
  isTrivial = \(s, i) -> if s == bit 0  && i == 1 then MinPlus 0 then Nothing,
  subproblems = \(s, i) ->
    if i == 1 then
      if s==bit n -1 then [(MinPlus $ dists!(j,1), (s, j)) | j <- [2..n]] 
      else []
    else
    let s' =  s .&. complement (bit (i-1)) in
      [(MinPlus $ dists!(i,j), (s', j)) | j <- [1..n], s' .&. bit (j-1) /= 0]
}

-- DoubleにはmaxBoundがないため、別途インスタンス化
instance Semiring MinPlus where
  t1@(MinPlusDouble v1) <+> t2@(MinPlusDouble v2) = MinPlusDouble (min v1 v2)
  t1@(MinPlusDouble v1) <.> t2@(MinPlusDouble v2)
    | t1 == zero = zero
    | t2 == zero = zero
    | otherwise = MinPlusDouble (v1 + v2)
  zero = MinPlusDouble 1e38
  one = MinPlusDouble 0

しゃくとり法+DP+剰余半環

ModIntはこちらのスクラップを参照。
しゃくとり法はこちらの記事を参照。

supplement :: (Array Int Int, Int, Int) -> DPProblem Int (ModInt 1000000007)
supplement (as,n,m) = DPProblem {
  start = n,
  getRange = (-1, n),
  isTrivial = \i-> case i of
    -1 -> Just (-1)
    0 -> Just 0
    1 -> Just 1
    _ -> Nothing,
  subproblems = \i -> [(2, i - 1), (-1, i-as ! i-1)]
}

Discussion