Haskell中的记忆？

Question

关于如何有效地解决Haskell中的以下函数的任何指针,对于大数 (n > 108)

f(n) = max(n, f(n/2) + f(n/3) + f(n/4))

我已经在Haskell中看到了用于解决斐波纳契数的例子,其中涉及计算(懒惰)所有斐波纳契数到所需的n.但在这种情况下,对于给定的n,我们只需要计算很少的中间结果.

谢谢

Answer 1

我们可以通过制作一个可以在亚线性时间内索引的结构来非常有效地完成这项工作.

但首先,

{-# LANGUAGE BangPatterns #-}

import Data.Function (fix)

让我们定义f,但让它使用'open recursion'而不是直接调用自己.

f :: (Int -> Int) -> Int -> Int
f mf 0 = 0
f mf n = max n $ mf (n `div` 2) +
                 mf (n `div` 3) +
                 mf (n `div` 4)

你可以f通过使用来获得一个unmemoizedfix f

这将允许您测试通过调用f执行您对小值的意义f,例如:fix f 123 = 144

我们可以通过定义来记住这个:

f_list :: [Int]
f_list = map (f faster_f) [0..]

faster_f :: Int -> Int
faster_f n = f_list !! n

这表现得非常好,并用记忆中间结果的东西取代了O(n ^ 3)时间.

但它仍然需要线性时间才能找到记忆的答案mf.这意味着结果如下:

*Main Data.List> faster_f 123801
248604

是可以容忍的,但结果并没有比这更好.我们可以做得更好!

首先,让我们定义一个无限树:

data Tree a = Tree (Tree a) a (Tree a)
instance Functor Tree where
    fmap f (Tree l m r) = Tree (fmap f l) (f m) (fmap f r)

然后我们将定义一种索引方式,因此我们可以找到索引n在O(log n)时间的节点:

index :: Tree a -> Int -> a
index (Tree _ m _) 0 = m
index (Tree l _ r) n = case (n - 1) `divMod` 2 of
    (q,0) -> index l q
    (q,1) -> index r q

...我们可能会发现一个充满自然数字的树很方便,所以我们不必乱用这些指数:

nats :: Tree Int
nats = go 0 1
    where
        go !n !s = Tree (go l s') n (go r s')
            where
                l = n + s
                r = l + s
                s' = s * 2

由于我们可以索引,您只需将树转换为列表:

toList :: Tree a -> [a]
toList as = map (index as) [0..]

到目前为止,您可以通过验证toList nats给您的工作来检查工作[0..]

现在,

f_tree :: Tree Int
f_tree = fmap (f fastest_f) nats

fastest_f :: Int -> Int
fastest_f = index f_tree

与上面的列表一样工作,但不是花费线性时间来查找每个节点,而是可以在对数时间内追逐它.

结果相当快:

*Main> fastest_f 12380192300
67652175206

*Main> fastest_f 12793129379123
120695231674999

事实上它是如此之快,以至于你可以通过以上方式进行替换Int,Integer并且几乎可以立即获得可笑的大答案

*Main> fastest_f' 1230891823091823018203123
93721573993600178112200489

*Main> fastest_f' 12308918230918230182031231231293810923
11097012733777002208302545289166620866358

Answer 2

爱德华的答案是如此美妙的宝石,我已经复制了它并提供了以开放递归形式记忆功能的组合memoList和memoTree组合器的实现.

{-# LANGUAGE BangPatterns #-}

import Data.Function (fix)

f :: (Integer -> Integer) -> Integer -> Integer
f mf 0 = 0
f mf n = max n $ mf (div n 2) +
                 mf (div n 3) +
                 mf (div n 4)


-- Memoizing using a list

-- The memoizing functionality depends on this being in eta reduced form!
memoList :: ((Integer -> Integer) -> Integer -> Integer) -> Integer -> Integer
memoList f = memoList_f
  where memoList_f = (memo !!) . fromInteger
        memo = map (f memoList_f) [0..]

faster_f :: Integer -> Integer
faster_f = memoList f


-- Memoizing using a tree

data Tree a = Tree (Tree a) a (Tree a)
instance Functor Tree where
    fmap f (Tree l m r) = Tree (fmap f l) (f m) (fmap f r)

index :: Tree a -> Integer -> a
index (Tree _ m _) 0 = m
index (Tree l _ r) n = case (n - 1) `divMod` 2 of
    (q,0) -> index l q
    (q,1) -> index r q

nats :: Tree Integer
nats = go 0 1
    where
        go !n !s = Tree (go l s') n (go r s')
            where
                l = n + s
                r = l + s
                s' = s * 2

toList :: Tree a -> [a]
toList as = map (index as) [0..]

-- The memoizing functionality depends on this being in eta reduced form!
memoTree :: ((Integer -> Integer) -> Integer -> Integer) -> Integer -> Integer
memoTree f = memoTree_f
  where memoTree_f = index memo
        memo = fmap (f memoTree_f) nats

fastest_f :: Integer -> Integer
fastest_f = memoTree f

Answer 3

不是最有效的方式,但记住:

f = 0 : [ g n | n <- [1..] ]
    where g n = max n $ f!!(n `div` 2) + f!!(n `div` 3) + f!!(n `div` 4)

在请求时f !! 144,检查是否f !! 143存在,但不计算其确切值.它仍然被设置为一些未知的计算结果.计算出的唯一精确值是所需的值.

所以最初,就计算了多少而言,程序一无所知.

f = .... 

当我们发出请求时f !! 12,它会开始进行一些模式匹配:

f = 0 : g 1 : g 2 : g 3 : g 4 : g 5 : g 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

现在开始计算了

f !! 12 = g 12 = max 12 $ f!!6 + f!!4 + f!!3

这递归地对f产生了另一个需求,所以我们计算

f !! 6 = g 6 = max 6 $ f !! 3 + f !! 2 + f !! 1
f !! 3 = g 3 = max 3 $ f !! 1 + f !! 1 + f !! 0
f !! 1 = g 1 = max 1 $ f !! 0 + f !! 0 + f !! 0
f !! 0 = 0

现在我们可以回流一些

f !! 1 = g 1 = max 1 $ 0 + 0 + 0 = 1

这意味着该程序现在知道:

f = 0 : 1 : g 2 : g 3 : g 4 : g 5 : g 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

继续涓涓细流:

f !! 3 = g 3 = max 3 $ 1 + 1 + 0 = 3

这意味着该程序现在知道:

f = 0 : 1 : g 2 : 3 : g 4 : g 5 : g 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

现在我们继续计算f!!6:

f !! 6 = g 6 = max 6 $ 3 + f !! 2 + 1
f !! 2 = g 2 = max 2 $ f !! 1 + f !! 0 + f !! 0 = max 2 $ 1 + 0 + 0 = 2
f !! 6 = g 6 = max 6 $ 3 + 2 + 1 = 6

这意味着该程序现在知道:

f = 0 : 1 : 2 : 3 : g 4 : g 5 : 6 : g 7 : g 8 : g 9 : g 10 : g 11 : g 12 : ...

现在我们继续计算f!!12:

f !! 12 = g 12 = max 12 $ 6 + f!!4 + 3
f !! 4 = g 4 = max 4 $ f !! 2 + f !! 1 + f !! 1 = max 4 $ 2 + 1 + 1 = 4
f !! 12 = g 12 = max 12 $ 6 + 4 + 3 = 13

这意味着该程序现在知道:

f = 0 : 1 : 2 : 3 : 4 : g 5 : 6 : g 7 : g 8 : g 9 : g 10 : g 11 : 13 : ...

所以计算是相当懒惰的.程序知道f !! 8存在的某些值,它等于g 8,但它不知道是什么g 8.

Answer 4

这是Edward Kmett出色答案的附录.

当我尝试他的代码时,定义nats和index看起来很神秘,所以我写了一个我发现更容易理解的替代版本.

我定义index并nats根据index'和nats'.

index' t n是在范围内定义的[1..].(回想一下,它index t是在范围内定义的[0..].)它可以通过将树n视为一串位来搜索树,然后反向读取这些位.如果该位是1,则采用右侧分支.如果该位是0,则采用左侧分支.它到达最后一位(必须是a 1)时停止.

index' (Tree l m r) 1 = m
index' (Tree l m r) n = case n `divMod` 2 of
                          (n', 0) -> index' l n'
                          (n', 1) -> index' r n'

正如nats定义的index那样index nats n == n始终为真,nats'定义为index'.

nats' = Tree l 1 r
  where
    l = fmap (\n -> n*2)     nats'
    r = fmap (\n -> n*2 + 1) nats'
    nats' = Tree l 1 r

现在,nats和index仅仅是nats'和index',但与1移值:

index t n = index' t (n+1)
nats = fmap (\n -> n-1) nats'

Answer 5

正如Edward Kmett的回答所述,为了加快速度,您需要缓存昂贵的计算并能够快速访问它们.

为了保持函数非monadic,构建无限懒惰树的解决方案,以及对其进行索引的适当方式(如先前帖子中所示)实现了该目标.如果放弃函数的非monadic性质,可以将Haskell中可用的标准关联容器与"状态"monad(如State或ST)结合使用.

虽然主要的缺点是你得到一个非monadic函数,你不必再自己索引结构,并且只能使用关联容器的标准实现.

为此,首先需要重新编写函数来接受任何类型的monad:

fm :: (Integral a, Monad m) => (a -> m a) -> a -> m a
fm _    0 = return 0
fm recf n = do
   recs <- mapM recf $ div n <$> [2, 3, 4]
   return $ max n (sum recs)

对于您的测试,您仍然可以使用Data.Function.fix定义一个不进行任何记忆的函数,尽管它有点冗长:

noMemoF :: (Integral n) => n -> n
noMemoF = runIdentity . fix fm

然后,您可以将State monad与Data.Map结合使用以加快速度:

import qualified Data.Map.Strict as MS

withMemoStMap :: (Integral n) => n -> n
withMemoStMap n = evalState (fm recF n) MS.empty
   where
      recF i = do
         v <- MS.lookup i <$> get
         case v of
            Just v' -> return v' 
            Nothing -> do
               v' <- fm recF i
               modify $ MS.insert i v'
               return v'

通过稍作更改,您可以调整代码以使用Data.HashMap:

import qualified Data.HashMap.Strict as HMS

withMemoStHMap :: (Integral n, Hashable n) => n -> n
withMemoStHMap n = evalState (fm recF n) HMS.empty
   where
      recF i = do
         v <- HMS.lookup i <$> get
         case v of
            Just v' -> return v' 
            Nothing -> do
               v' <- fm recF i
               modify $ HMS.insert i v'
               return v'

您也可以尝试将可变数据结构(如Data.HashTable)与ST monad结合使用,而不是持久数据结构:

import qualified Data.HashTable.ST.Linear as MHM

withMemoMutMap :: (Integral n, Hashable n) => n -> n
withMemoMutMap n = runST $
   do ht <- MHM.new
      recF ht n
   where
      recF ht i = do
         k <- MHM.lookup ht i
         case k of
            Just k' -> return k'
            Nothing -> do 
               k' <- fm (recF ht) i
               MHM.insert ht i k'
               return k'

与没有任何memoization的实现相比,任何这些实现​​都允许您在巨大的输入下以微秒为单位获得结果,而不必等待几秒钟.

使用Criterion作为基准,我可以观察到Data.HashMap的实现实际上比定时非常相似的Data.Map和Data.HashTable略好(大约20%).

我发现基准测试的结果有点令人惊讶.我最初的感觉是HashTable的性能优于HashMap,因为它是可变的.在最后一个实现中可能隐藏了一些性能缺陷.