¿Cómo se escriben algoritmos de programación dinámica eficientes en Haskell?

Question

He estado jugando con programación dinámica en Haskell. Prácticamente cada tutorial que he visto sobre el tema ofrece el mismo algoritmo muy elegante basado en la memorización y la pereza del tipo Array. Inspirado por esos ejemplos, escribí el siguiente algoritmo como una prueba:

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

Mi único problema es la eficiencia. Incluso usando -O2 de GHC, este programa demora 1.6 segundos para calcular pascal 1000, que es aproximadamente 160 veces más lento que un programa C++ no optimizado equivalente. Y la brecha solo se amplía con entradas más grandes.

Parece que he intentado todas las permutaciones posibles del código anterior, junto con las alternativas sugeridas, como la biblioteca de datacombinators, y todas tuvieron el mismo o peor rendimiento. Lo único que no he intentado es la ST Monad, que estoy seguro se podría hacer para ejecutar el programa solo más ligero más lento que la versión C. Pero me gustaría escribirlo en el idiomático Haskell, y no entiendo por qué la versión idiomática es tan ineficiente. Tengo dos preguntas:

1. ¿Por qué el código anterior es tan ineficiente? Parece una iteración directa a través de una matriz, con una operación aritmética en cada entrada. Claramente, Haskell está haciendo algo detrás de escena que no entiendo.

2. Hay una manera de hacerlo mucho más eficiente (a lo sumo 10-15 veces el tiempo de ejecución de un programa C) sin sacrificar su formulación recursiva sin estado (en comparación con una implementación que utiliza matrices mutables en el ST Monada)?

Muchas gracias.

Editar: El módulo de matriz utilizado es el estándar Data.Array

Answer 1

Bueno, el algoritmo podría diseñarse un poco mejor . Usando el paquete vector y ser inteligentes acerca de mantener sólo una fila en la memoria a la vez, podemos conseguir algo que es idiomática de una manera diferente:

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

Esto optimiza abajo muy fuertemente, sobre todo porque el paquete vector incluye algunos trucos bastante ingeniosos para optimizar de forma transparente las operaciones de matriz escritas en un estilo idiomático.

Answer 2

1 ¿Por qué es el código anterior tan ineficiente? Parece una iteración directa a través de una matriz, con una operación aritmética en cada entrada. Claramente, Haskell está haciendo algo detrás de escena que no entiendo.

El problema es que el código escribe thunks en la matriz. Luego, cuando se lee la entrada (n,n), la evaluación de los thunks salta de nuevo por toda la matriz, recurriendo hasta que finalmente se encuentra un valor que no necesita más recursión. Eso causa mucha asignación e ineficiencia innecesarias.

El código C++ no tiene ese problema, los valores se escriben y se leen directamente sin requerir una evaluación adicional. Como ocurriría con un STUArray. ¿Tiene

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 

realmente se ven tan mal?

2 ¿Hay una manera de hacer que sea mucho más eficiente (en la mayoría de 10-15 veces el tiempo de ejecución de un programa en C) sin sacrificar su formulación sin estado, recursiva (vis-a-vis una implementación utilizando matrices en mutables la ST Monad)?

No sé de uno. Pero puede haber.

Adición:

Una vez que uno utiliza STUArray s o sin embalaje Vector s, todavía hay una diferencia significativa a la implementación en C equivalente. La razón es que gcc reemplaza el % por una combinación de multiplicaciones, desplazamientos y sustracciones (incluso sin optimizaciones), ya que el módulo es conocido. Hacer lo mismo con la mano en Haskell (GHC ya no [todavía] hacer eso),

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

consigue las versiones Haskell a la par con C.

Answer 3

El truco es pensar cómo escribir todo el maldito algoritmo de una vez, y luego usar vectores sin caja como su tipo de datos de respaldo. Por ejemplo, la siguiente funciona con cerca de 20 veces más rápido en mi máquina de su código:

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)) 

Entonces escribí dos main funciones que llamó a la suya y la mía con un argumento de 4000; estos funcionaron en 10.42s y 0.54s respectivamente. Por supuesto, como estoy seguro de que sabe, ambos se vuelen fuera del agua (0.00s) por la versión que utiliza un algoritmo mejor:

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