Memoization w Haskell?

Odpowiedź Edwarda jest tak cudownym klejnotem, że zduplikowałem ją i dostarczyłem implementacji memoList i memoTree kombinatorów, które zapamiętują funkcję w formie otwartej rekurencyjnej.

{-# 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

Nie najskuteczniejszy sposób, ale robi memoize:

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

Podczas żądania f !! 144 sprawdza się, czy f !! 143 istnieje, ale jego dokładna wartość nie jest obliczana. Nadal jest to jakiś nieznany wynik obliczeń. Tylko dokładne obliczone wartości są potrzebne.

Więc początkowo, jeśli chodzi o to, ile zostało obliczone, program nic nie wie.

f = .... 

Kiedy wykonamy żądanie f !! 12, zaczyna ono dopasowywać wzór:

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 : ...

Teraz to rozpoczyna obliczanie

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

To rekurencyjnie powoduje kolejne zapotrzebowanie na f, więc obliczamy

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

Teraz możemy zrobić backup

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

Co oznacza, że program już wie:

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 : ...

Ciąg dalszy:

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

Co oznacza, że program już wie:

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

Teraz kontynuujemy nasze obliczenia 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

Co oznacza, że program już wie:

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

Teraz kontynuujemy nasze obliczenia 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

Co oznacza, że program już wie:

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

Więc obliczenia są wykonywane dość leniwie. Program wie, że istnieje jakaś wartość f !! 8, która jest równa g 8, ale nie ma pojęcia, czym jest g 8.

Jak stwierdzono w odpowiedzi Edwarda Kmetta, aby przyspieszyć wszystko, musisz buforować kosztowne obliczenia i mieć do nich szybki dostęp.

Aby funkcja nie była monadyczna, rozwiązanie budowania nieskończonego leniwego drzewa z odpowiednim sposobem indeksowania go (jak pokazano w poprzednich postach) spełnia ten cel. Jeśli zrezygnujesz z niemonadowego charakteru funkcji, możesz użyć standardowych kontenerów asocjacyjnych dostępnych w Haskell w połączeniu z monadami "stanowymi" (takimi jak stan lub ST).

Podczas gdy główną wadą jest to, że otrzymujesz funkcję niemonadyczną, nie musisz już sam indeksować struktury i możesz po prostu korzystać ze standardowych implementacji kontenerów asocjacyjnych.

Aby to zrobić, musisz najpierw napisać ponownie funkcję, aby zaakceptować dowolny rodzaj monady:

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)

Dla Twoich testów, nadal możesz zdefiniować funkcję, która nie robi memoizacji za pomocą danych.Funkcja.fix, chociaż jest nieco bardziej wyrazisty:

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

Możesz wtedy użyć State monad w połączenie z danymi.Mapa przyspieszająca:

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'

Dzięki niewielkim zmianom można dostosować kod do pracy z danymi.HashMap zamiast:

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'

Zamiast trwałych struktur danych, możesz również wypróbować zmienne struktury danych (takie jak dane.HashTable) w połączeniu ze 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'

W porównaniu z implementacją bez żadnej memoizacji, każda z tych implementacji pozwala na uzyskanie wyników w mikro-sekundach, zamiast czekać kilka sekund.

Używając kryterium jako benchmark, mogłem zauważyć, że implementacja z danymi.HashMap faktycznie wypadł nieco lepiej (około 20%) niż dane.Mapa i dane.HashTable, dla których czasy były bardzo podobne.

Wyniki benchmarka uznałem za nieco zaskakujące. Początkowo miałem wrażenie, że HashTable przewyższy implementację HashMap, ponieważ jest zmienna. W tej ostatniej implementacji może być ukryty pewien defekt wydajności.

To jest dodatek do doskonałej odpowiedzi Edwarda Kmetta.

Kiedy próbowałem jego kodu, definicje nats i index wydawały się dość tajemnicze, więc napisałem alternatywną wersję, która okazała się łatwiejsza do zrozumienia.

Definiuję index i nats w kategoriach index' i nats'.

index' t n jest zdefiniowany w zakresie [1..]. (Przypomnijmy, że {[11] } jest zdefiniowany w zakresie [0..].) Działa przeszukiwanie drzewa traktując n jako ciąg bitów i czytając bity w odwrotnej kolejności. Jeśli bit jest 1, to pobiera gałąź prawej ręki. Jeśli bit jest 0, przyjmuje gałąź lewej ręki. Zatrzymuje się, gdy osiągnie ostatni bit (który musi być 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'

Tak jak nats jest zdefiniowane dla index tak, że index nats n == n jest zawsze prawdziwe, nats' jest zdefiniowane dla 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

Teraz nats i index są po prostu nats' i index', ale z wartościami przesuniętymi o 1:

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

Kilka lat później, spojrzałem na to i zdałem sobie sprawę, że istnieje prosty sposób na zapamiętanie tego w czasie liniowym za pomocą zipWith i funkcji pomocniczej:

dilate :: Int -> [x] -> [x]
dilate n xs = replicate n =<< xs

dilate posiada właściwość handy dilate n xs !! i == xs !! div i n.

Więc, przypuśćmy, że otrzymamy f (0), upraszcza to obliczenie do

fs = f0 : zipWith max [1..] (tail $ fs#/2 .+. fs#/3 .+. fs#/4)
  where (.+.) = zipWith (+)
        infixl 6 .+.
        (#/) = flip dilate
        infixl 7 #/

Wyglądający bardzo podobnie do naszego oryginalnego opisu problemu i dający liniowe rozwiązanie (sum $ take n fs zajmie O (n)).

Kolejny dodatek do odpowiedzi Edwarda Kmetta: samodzielny przykład:

data NatTrie v = NatTrie (NatTrie v) v (NatTrie v)

memo1 arg_to_index index_to_arg f = (\n -> index nats (arg_to_index n))
  where nats = go 0 1
        go i s = NatTrie (go (i+s) s') (f (index_to_arg i)) (go (i+s') s')
          where s' = 2*s
        index (NatTrie l v r) i
          | i <  0    = f (index_to_arg i)
          | i == 0    = v
          | otherwise = case (i-1) `divMod` 2 of
             (i',0) -> index l i'
             (i',1) -> index r i'

memoNat = memo1 id id 

Użyj go w następujący sposób, aby zapamiętać funkcję z pojedynczą liczbą całkowitą arg (np.]}

fib = memoNat f
  where f 0 = 0
        f 1 = 1
        f n = fib (n-1) + fib (n-2)

Tylko wartości dla nieujemnych argumentów będą buforowane.

Aby również buforować wartości ujemnych argumentów, użyj memoInt, zdefiniowanego w następujący sposób:

memoInt = memo1 arg_to_index index_to_arg
  where arg_to_index n
         | n < 0     = -2*n
         | otherwise =  2*n + 1
        index_to_arg i = case i `divMod` 2 of
           (n,0) -> -n
           (n,1) ->  n

Do wartości bufora dla funkcji z dwoma argumentami całkowitymi użyj memoIntInt, zdefiniowanego w następujący sposób:

memoIntInt f = memoInt (\n -> memoInt (f n))

Rozwiązanie bez indeksowania, a nie oparte na Edwardzie KMETCIE.]}

Uwzględniam wspólne podzbiory do wspólnego rodzica (f(n/4) jest dzielone pomiędzy f(n/2) i f(n/4), A f(n/6) jest dzielone między f(2) i f(3)). Zapisując je jako pojedynczą zmienną w rodzicu, obliczenie poddrzewa jest wykonywane raz.

data Tree a =
  Node {datum :: a, child2 :: Tree a, child3 :: Tree a}

f :: Int -> Int
f n = datum root
  where root = f' n Nothing Nothing


-- Pass in the arg
  -- and this node's lifted children (if any).
f' :: Integral a => a -> Maybe (Tree a) -> Maybe (Tree a)-> a
f' 0 _ _ = leaf
    where leaf = Node 0 leaf leaf
f' n m2 m3 = Node d c2 c3
  where
    d = if n < 12 then n
            else max n (d2 + d3 + d4)
    [n2,n3,n4,n6] = map (n `div`) [2,3,4,6]
    [d2,d3,d4,d6] = map datum [c2,c3,c4,c6]
    c2 = case m2 of    -- Check for a passed-in subtree before recursing.
      Just c2' -> c2'
      Nothing -> f' n2 Nothing (Just c6)
    c3 = case m3 of
      Just c3' -> c3'
      Nothing -> f' n3 (Just c6) Nothing
    c4 = child2 c2
    c6 = f' n6 Nothing Nothing

    main =
      print (f 123801)
      -- Should print 248604.

Kod nie łatwo rozszerza się na ogólną funkcję memoizacji (przynajmniej nie wiedziałbym, jak to zrobić), a naprawdę musisz pomyśleć, jak podproblemy nakładają się na siebie, ale strategia powinna działać dla ogólnych wielu nie-całkowitych parametrów. (Wymyśliłem to dla dwóch parametrów strunowych.)

Notatka jest odrzucana po każdym obliczeniu. (Znowu myślałem o dwóch parametrach strunowych.)

Nie wiem, czy to jest bardziej skuteczne niż inne odpowiedzi. Każde wyszukiwanie jest technicznie tylko jednym lub dwoma krokami ("spójrz na swoje dziecko lub dziecko dziecka"), ale może być dużo dodatkowego wykorzystania pamięci.

Edit: To rozwiązanie jeszcze się nie zgadza. Podział jest niekompletny.

Edit: teraz powinno to być dzielenie się dziećmi, ale zdałem sobie sprawę, że ten problem ma wiele nietrywialnego dzielenia się: n/2/2/2 i n/3/3 mogą być takie same. Problem nie pasuje do mojej strategii.