Quadratic complexity though use of STArrays

Dear haskell-cafe users, I am constructing a shuffle function: given an StdGen and a list, return the list permuted, with all permutations of equal probability. There is the simlpe recursive definition: generate a number from 1 to length list, take this element out from the list, call the function recursively on the remaining list and then cons the element on the shuffled list. A more imperative approach is to make the list an array, and traverse the array in reverse, swapping the iterated element with an arbitrary element less than or equal to the iterator. These functions are implemented as shuffleRec and shuffleArr, respectively. What complexity does these functions have? I argue that the shuffleArr function should be O(n), since it only contains one loop of n, where each loop does actions that are O(1): generating a random number and swapping two elements in an array. I argue that the shuffleRec function should be O(n^2), since for each call, it creates a new list in O(n), with the drop and take calls, and calls itself recursively. This yields O(n^2). However, they both have the same runnig time (roughly), and through looking at the plot it _very_ much looks quadratic. I am compiling with GHC and I guess there is something in the lazy semantics that confuses me about the complexities, and maybe I have misunderstood how STArrays work. Any pointers to what's going in is greatly appreciated! Best regards, Dan Rosén, Sweden Here is the code: module Main where import Control.Monad import Control.Monad.ST import Data.Array.ST import Data.STRef import System.Random import Time import CPUTime shuffleArr :: StdGen -> [a] -> [a] shuffleArr g list = runST $ do let n = length list gref <- newSTRef g arr <- listToArray list forM_ [n,n-1..2] $ \p -> do m <- rand (1,p) gref swap arr m p getElems arr where rand range gref = do g <- readSTRef gref let (v,g') = randomR range g writeSTRef gref g' return v swap a n m = do [n',m'] <- mapM (readArray a) [n,m] mapM (uncurry $ writeArray a) [(m,n'),(n,m')] listToArray :: [a] -> ST s (STArray s Int a) listToArray list = let n = length list in newListArray (1,n) list shuffleRec :: StdGen -> [a] -> [a] shuffleRec g list = x:shuffleArr g' xs where (n,g') = randomR (0,length list-1) g (x:xs') = drop n list xs = take n list ++ xs' -- A somewhat lame attempt to derive the complexities through testing, -- prints the times for the different functions in a table main :: IO () main = do let times = take 30 $ iterate (+30000) 10000 answers <- mapM test times sequence_ [ putStrLn $ concatMap ((++ "\t"). show) [toInteger t,arr,rec] | (t,(arr,rec)) <- zip times answers ] -- Perform a test of size n, and return the cycles it took for the different -- algorithms in a pair. Evaluation is enforced by seq on length of the list. test :: Int -> IO (Integer,Integer) test n = do let list = [1..n] [g1,g2] <- replicateM 2 newStdGen length list `seq` do s <- doTime ("shuffleArr " ++ show n) $ (length $ shuffleArr g1 list) `seq` return () s' <- doTime ("shuffleRec " ++ show n) $ (length $ shuffleRec g2 list) `seq` return () return (s,s') -- This is taken from GenUtil from the JHC creator's homepage doTime :: String -> IO a -> IO Integer doTime str action = do start <- getCPUTime x <- action end <- getCPUTime let time = (end - start) `div` 1000000 -- `div` cpuTimePrecision -- putStrLn $ "Timing: " ++ str ++ " " ++ show time return time

Am Dienstag 22 September 2009 21:31:08 schrieb Dan Rosén:
Dear haskell-cafe users,
I am constructing a shuffle function: given an StdGen and a list, return the list permuted, with all permutations of equal probability.
There is the simlpe recursive definition: generate a number from 1 to length list, take this element out from the list, call the function recursively on the remaining list and then cons the element on the shuffled list.
A more imperative approach is to make the list an array, and traverse the array in reverse, swapping the iterated element with an arbitrary element less than or equal to the iterator.
These functions are implemented as shuffleRec and shuffleArr, respectively.
What complexity does these functions have?
I argue that the shuffleArr function should be O(n), since it only contains one loop of n, where each loop does actions that are O(1): generating a random number and swapping two elements in an array.
I argue that the shuffleRec function should be O(n^2), since for each call, it creates a new list in O(n), with the drop and take calls, and calls itself recursively. This yields O(n^2).
However, they both have the same runnig time (roughly), and through looking at the plot it _very_ much looks quadratic.
Regarding
shuffleRec :: StdGen -> [a] -> [a] shuffleRec g list = x:shuffleArr g' xs where (n,g') = randomR (0,length list-1) g (x:xs') = drop n list xs = take n list ++ xs'
it's not surprising they take more or less the same time. Make it shuffleRec g list = x:shuffleRec g' xs and prepare to kill the process pretty soon. That doesn't explain the quadratic behaviour of shuffleArr, though. I suspect it's laziness, things aren't actually done until the result is finally demanded, but I would have to take a closer look to really find out.
I am compiling with GHC and I guess there is something in the lazy semantics that confuses me about the complexities, and maybe I have misunderstood how STArrays work.
Any pointers to what's going in is greatly appreciated!
Best regards, Dan Rosén, Sweden
Here is the code:

Am Dienstag 22 September 2009 22:31:48 schrieb Daniel Fischer:
That doesn't explain the quadratic behaviour of shuffleArr, though. I suspect it's laziness, things aren't actually done until the result is finally demanded, but I would have to take a closer look to really find out.
Yep. Strictifying things gives the expected linear behaviour. In swap a n m = do [n',m'] <- mapM (readArray a) [n,m] mapM (uncurry $ writeArray a) [(m,n'),(n,m')] you don't actually read and write values, but ever longer thunks. Changing swap to swap a m n vm <- readArray a m vn <- readArray a n writeArray a n $! vm writeArray a m $! vn is all you need to do (as long as your values have a simple type).

Thanks Daniel and Tobias, and sorry for posting something with an obvious miss in it (the function call) It works better with strict evaulation. best regards, Dan

Hi Dan!
You might want to change the following:
shuffleRec :: StdGen -> [a] -> [a]
shuffleRec g list = x:shuffleArr g' xs
where
(n,g') = randomR (0,length list-1) g
(x:xs') = drop n list
xs = take n list ++ xs'
into the following:
shuffleRec :: StdGen -> [a] -> [a]
shuffleRec g list = x:shuffleRec g' xs
where
(n,g') = randomR (0,length list-1) g
(x:xs') = drop n list
xs = take n list ++ xs'
Since shuffleRec just called shuffleArr, one would expect them
to run in approximately the same time :-)
//Tobias
2009/9/22 Dan Rosén
Dear haskell-cafe users,
I am constructing a shuffle function: given an StdGen and a list, return the list permuted, with all permutations of equal probability.
There is the simlpe recursive definition: generate a number from 1 to length list, take this element out from the list, call the function recursively on the remaining list and then cons the element on the shuffled list.
A more imperative approach is to make the list an array, and traverse the array in reverse, swapping the iterated element with an arbitrary element less than or equal to the iterator.
These functions are implemented as shuffleRec and shuffleArr, respectively.
What complexity does these functions have?
I argue that the shuffleArr function should be O(n), since it only contains one loop of n, where each loop does actions that are O(1): generating a random number and swapping two elements in an array.
I argue that the shuffleRec function should be O(n^2), since for each call, it creates a new list in O(n), with the drop and take calls, and calls itself recursively. This yields O(n^2).
However, they both have the same runnig time (roughly), and through looking at the plot it _very_ much looks quadratic.
I am compiling with GHC and I guess there is something in the lazy semantics that confuses me about the complexities, and maybe I have misunderstood how STArrays work.
Any pointers to what's going in is greatly appreciated!
Best regards, Dan Rosén, Sweden
Here is the code:
module Main where
import Control.Monad import Control.Monad.ST import Data.Array.ST import Data.STRef import System.Random
import Time import CPUTime
shuffleArr :: StdGen -> [a] -> [a] shuffleArr g list = runST $ do let n = length list gref <- newSTRef g arr <- listToArray list forM_ [n,n-1..2] $ \p -> do m <- rand (1,p) gref swap arr m p getElems arr where rand range gref = do g <- readSTRef gref let (v,g') = randomR range g writeSTRef gref g' return v
swap a n m = do [n',m'] <- mapM (readArray a) [n,m] mapM (uncurry $ writeArray a) [(m,n'),(n,m')]
listToArray :: [a] -> ST s (STArray s Int a) listToArray list = let n = length list in newListArray (1,n) list
shuffleRec :: StdGen -> [a] -> [a] shuffleRec g list = x:shuffleArr g' xs where (n,g') = randomR (0,length list-1) g (x:xs') = drop n list xs = take n list ++ xs'
-- A somewhat lame attempt to derive the complexities through testing, -- prints the times for the different functions in a table main :: IO () main = do let times = take 30 $ iterate (+30000) 10000 answers <- mapM test times sequence_ [ putStrLn $ concatMap ((++ "\t"). show) [toInteger t,arr,rec] | (t,(arr,rec)) <- zip times answers ]
-- Perform a test of size n, and return the cycles it took for the different -- algorithms in a pair. Evaluation is enforced by seq on length of the list. test :: Int -> IO (Integer,Integer) test n = do let list = [1..n] [g1,g2] <- replicateM 2 newStdGen length list `seq` do s <- doTime ("shuffleArr " ++ show n) $ (length $ shuffleArr g1 list) `seq` return () s' <- doTime ("shuffleRec " ++ show n) $ (length $ shuffleRec g2 list) `seq` return () return (s,s')
-- This is taken from GenUtil from the JHC creator's homepage doTime :: String -> IO a -> IO Integer doTime str action = do start <- getCPUTime x <- action end <- getCPUTime let time = (end - start) `div` 1000000 -- `div` cpuTimePrecision -- putStrLn $ "Timing: " ++ str ++ " " ++ show time return time _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
-- Tobias Olausson tobsan@gmail.com

Dan Rosén wrote:
What complexity does these functions have?
I argue that the shuffleArr function should be O(n), since it only contains one loop of n, where each loop does actions that are O(1): generating a random number and swapping two elements in an array.
However, they both have the same runnig time (roughly), and through looking at the plot it _very_ much looks quadratic.
Right. In the case of shuffleArr, it's the garbage collection time that explodes. I think that the reason is that the garbage collector has to scan the whole array (which lives in the old generation) on each garbage collection, while the young generation is kept very small, making GCs very frequent. Indeed the program behaves much better with just one GC generation:
./a.out +RTS -G1 (1000000, 1988698, 1760732, 1.98870,1.76073) (1500000, 2991546, 2683592, 1.99436,1.78906) (2000000, 4018388, 3611451, 2.00919,1.80573)
Output format: (n, time for your original shuffleArray (t1), time for Daniel's stricter version (t2), t1/n, t2/n) While with the default settings, it looks much worse:
./a.out (1000000, 7575850, 7958790, 7.57585, 7.95879) (1500000,17119397,19449044,11.41293,12.96602) (2000000,30687335,35125661,15.34367,17.56283)
You can also increase the initial heap size,
./a.out +RTS -H250M (1000000, 1131828, 901863, 1.13183, 0.90186) (1500000, 1726737, 1427783, 1.15116, 0.95186) (2000000, 2324647, 1935705, 1.16232, 0.96785)
Bertram
participants (5)
-
Bertram Felgenhauer
-
Dan Rosén
-
Dan Rosén
-
Daniel Fischer
-
Tobias Olausson