Come si scrivono algoritmi di programmazione dinamica efficienti in Haskell?

Question

Ho giocato con la programmazione dinamica in Haskell. Praticamente ogni tutorial che ho visto sull'argomento fornisce lo stesso algoritmo molto elegante basato sulla memoizzazione e la pigrizia del tipo di matrice. Ispirato da questi esempi, ho scritto il seguente algoritmo come un test:

-- pascal n returns the nth entry on the main diagonal of pascal's triangle 
-- (mod a million for efficiency) 
pascal :: Int -> Int 
pascal n = p ! (n,n) where 
      p = listArray ((0,0),(n,n)) [f (i,j) | i <- [0 .. n], j <- [0 .. n]] 

      f :: (Int,Int) -> Int 
      f (_,0) = 1 
      f (0,_) = 1 
      f (i,j) = (p ! (i, j-1) + p ! (i-1, j)) `mod` 1000000 

Il mio unico problema è l'efficienza. Anche usando GHC's -O2, questo programma impiega 1,6 secondi per calcolare pascal 1000, che è circa 160 volte più lento di un equivalente programma C++ non ottimizzato. E il divario si allarga solo con input più grandi.

Sembra che abbia provato ogni possibile permutazione del codice precedente, insieme a alternative suggerite come la libreria di dati-memocombinatori, e tutte hanno avuto le stesse o peggiori prestazioni. L'unica cosa che non ho provato è la ST Monad, che sono sicuro che potrebbe essere fatta per eseguire il programma solo più lentamente più lento della versione C. Ma mi piacerebbe davvero scriverlo nell'idiomatico Haskell, e non capisco perché la versione idiomatica sia così inefficiente. Ho due domande:

1. Perché il codice di cui sopra è così inefficiente? Sembra una semplice iterazione attraverso una matrice, con un'operazione aritmetica ad ogni voce. Chiaramente Haskell sta facendo qualcosa dietro le quinte che non capisco.

2. C'è un modo per renderlo molto più efficiente (al massimo 10-15 volte il tempo di esecuzione di un programma C) senza sacrificare la sua formulazione ricorsiva senza stato (di fronte a un'implementazione che utilizza matrici mutevoli nella ST Monade)?

Grazie mille.

Edit: Il modulo matrice utilizzata è lo standard Data.Array

Answer 1

Ebbene, l'algoritmo potrebbe essere progettato un po 'meglio . Usando il pacchetto vector e di essere intelligente su come mantenere solo una riga in memoria alla volta, possiamo ottenere qualcosa che è idiomatica in un modo diverso:

{-# LANGUAGE BangPatterns #-} 
import Data.Vector.Unboxed 
import Prelude hiding (replicate, tail, scanl) 

pascal :: Int -> Int 
pascal !n = go 1 ((replicate (n+1) 1) :: Vector Int) where 
    go !i !prevRow 
    | i <= n = go (i+1) (scanl f 1 (tail prevRow)) 
    | otherwise = prevRow ! n 
    f x y = (x + y) `rem` 1000000 

Ciò consente di ottimizzare le molto strettamente, soprattutto perché il pacchetto vector include alcuni trucchi piuttosto ingegnosi per ottimizzare in modo trasparente le operazioni dell'array scritte in uno stile idiomatico.

Answer 2

1 Perché il codice di cui sopra in modo inefficiente? Sembra una semplice iterazione attraverso una matrice, con un'operazione aritmetica ad ogni voce. Chiaramente Haskell sta facendo qualcosa dietro le quinte che non capisco.

Il problema è che il codice scrive i thunk nell'array. Quindi, quando viene letto l'articolo (n,n), la valutazione dei thunks salta di nuovo su tutto l'array, ricorrendo fino a quando non viene trovato un valore che non richiede ulteriore ricorsione. Ciò causa un sacco di allocazioni e inefficienze non necessarie.

Il codice C++ non presenta questo problema, i valori vengono scritti e letti direttamente senza richiedere ulteriori valutazioni. Come accadrebbe con uno STUArray.

p = runSTUArray $ do 
    arr <- newArray ((0,0),(n,n)) 1 
    forM_ [1 .. n] $ \i -> 
     forM_ [1 .. n] $ \j -> do 
      a <- readArray arr (i,j-1) 
      b <- readArray arr (i-1,j) 
      writeArray arr (i,j) $! (a+b) `rem` 1000000 
    return arr 

sembra davvero così male?

2 C'è un modo per rendere molto più efficiente (al massimo 10-15 volte il tempo di esecuzione di un programma C) senza sacrificare la sua stateless, formulazione ricorsiva (vis-a-vis un'implementazione utilizzando matrici mutabili in la Monad ST)?

Non so di uno. Ma potrebbe esserci.

Addendum:

Una volta che si usa STUArray s o disimballati Vector s, c'è ancora una differenza significativa per l'attuazione equivalente C. Il motivo è che gcc sostituisce lo % con una combinazione di moltiplicazioni, spostamenti e sottrazioni (anche senza ottimizzazioni), poiché il modulo è noto. Fare la stessa mano in Haskell (dal momento che non si GHC [ancora] farlo),

-- fast modulo 1000000 
-- for nonnegative Ints < 2^31 
-- requires 64-bit Ints 
fastMod :: Int -> Int 
fastMod n = n - 1000000*((n*1125899907) `shiftR` 50) 

ottiene le versioni Haskell alla pari con C.

Answer 3

Il trucco è pensare a come scrivere l'intero dannato algoritmo in una sola volta e quindi utilizzare i vettori non in scatola come tipo di dati di supporto. Ad esempio, il seguente gestisce circa 20 volte più veloce sulla mia macchina di tuo codice:

import qualified Data.Vector.Unboxed as V 

combine :: Int -> Int -> Int 
combine x y = (x+y) `mod` 1000000 

pascal n = V.last $ go n where 
    go 0 = V.replicate (n+1) 1 
    go m = V.scanl1 combine (go (m-1)) 

poi ho scritto due main funzioni che chiamò tua e la mia con un argomento di 4000; queste corse rispettivamente in 10.42s e 0.54s. Naturalmente, come io sono sicuro che voi sapete, entrambi vengono soffiate fuori dall'acqua (0.00s) per la versione che utilizza un algoritmo migliore:

pascal' :: Integer -> Integer 
pascal :: Int -> Int 
pascal' n = product [n+1..n*2] `div` product [2..n] 
pascal = fromIntegral . (`mod` 1000000) . pascal' . fromIntegral