Как написать эффективные алгоритмы динамического программирования в Haskell?

Я играл с динамическим программированием в Haskell. Практически каждый учебник, который я видел по этому вопросу, дает тот же, очень элегантный алгоритм, основанный на memoization и лень типа Array. Вдохновленный этими примерами, я написал в качестве теста следующий алгоритм:

-- pascal n returns the nth entry on the main diagonal of pascal 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

Моя единственная проблема - эффективность. Даже используя GHC-O2, эта программа занимает 1,6 секунды для вычисления pascal 1000, что примерно в 160 раз медленнее, чем эквивалентная неоптимизированная программа на С++. И разрыв только расширяется с большими входами.

Кажется, я пробовал все возможные перестановки вышеприведенного кода вместе с предложенными альтернативами, такими как библиотека данных-воспоминаний, и все они имели одинаковую или худшую производительность. Единственное, что я не пробовал, это ST Monad, и я уверен, что можно было запустить программу только медленнее, чем версия C. Но я действительно хотел бы написать его в идиоматическом Haskell, и я не понимаю, почему идиоматическая версия настолько неэффективна. У меня есть два вопроса:

Почему приведенный выше код настолько неэффективен? Это кажется простой итерацией через матрицу с арифметической операцией в каждой записи. Ясно, что Хаскелл делает что-то за кулисами, которые я не понимаю.
Есть ли способ сделать его намного более эффективным (не более 10-15 раз времени выполнения программы C), не жертвуя своей безрепортивной рекурсивной формулировкой (по сравнению с реализацией с использованием изменяемых массивов в ST Monad)?

Большое спасибо.

Изменить: используемый модуль массива - это стандартный Data.Array

Ответ 1

Ну, алгоритм может быть разработан немного лучше. Используя пакет vector и умеющий хранить только одну строку в памяти за раз, мы можем получить что-то идиоматическое по-другому:

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

Это оптимизируется очень сильно, особенно потому, что пакет vector содержит некоторые довольно хитроумные трюки для прозрачной оптимизации операций массива, написанных в идиоматическом стиле.

Ответ 2

1 Почему приведенный выше код настолько неэффективен? Это кажется простой итерацией через матрицу с арифметической операцией в каждой записи. Ясно, что Хаскелл делает что-то за кулисами, которые я не понимаю.

Проблема заключается в том, что код записывает thunks в массив. Затем, когда запись (n,n) считывается, оценка thunks снова перескакивает по всему массиву, повторяясь до тех пор, пока не будет найдено значение, не требующее дальнейшей рекурсии. Это вызывает много ненужного распределения и неэффективности.

Код С++ не имеет этой проблемы, значения записываются и читаются напрямую, не требуя дальнейшей оценки. Как это было бы с 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

действительно выглядит так плохо?

2 Есть ли способ сделать его намного более эффективным (не более 10-15 раз времени выполнения программы на C), не жертвуя своей безрепортивной рекурсивной формулировкой (по сравнению с реализацией с использованием изменяемых массивов в ST Monad )?

Я не знаю одного. Но может быть.

Приложение:

Как только один использует STUArray или unboxed Vector s, все еще значительная разница с эквивалентной реализацией C. Причина в том, что gcc заменяет % комбинацией умножений, сдвигов и вычитаний (даже без оптимизаций), поскольку модуль известен. Выполняя то же самое вручную в Haskell (поскольку GHC еще не делает этого),

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

получает версии Haskell на уровне с C.

Ответ 3

Трюк состоит в том, чтобы подумать о том, как написать весь проклятый алгоритм сразу, а затем использовать unboxed векторы в качестве вашего типа данных поддержки. Например, следующие программы работают примерно на 20 раз быстрее, чем ваш код:

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

Затем я написал две функции main, которые вызвали ваш и мой с аргументом 4000; они выполнялись в 10.42s и 0.54s соответственно. Конечно, как я уверен, вы знаете, они оба вырываются из воды (0.00s) версией, которая использует лучший алгоритм:

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